10 Most granular materials encountered in nature and industry lie either in the quasi-static regime or the intermediate dense ﬂow regime. Debris materials are a typical granular material with viscous interstitial ﬂuid, and shows solid-like behaviors before failure and ﬂuid-like behaviors after failure. Based on Bagnold’s pioneering work on granular-ﬂuid ﬂows, we propose a framework for constitutive model development, which has an ad- ditive form. Based on this framework, a uniﬁed constitutive model for granular-ﬂuid material in the quasi-static and dense ﬂow regimes is developed. The main intergranu- lar interactions and granular-ﬂuid interactions controlling the mechanical behaviors are taken into account using the Mohr-Coulomb model and a Bagnold-type relation. Dry granular ﬂows in three simple conﬁgurations, i.e., plain shear, vertical chute ﬂow and ﬂow on an inclined plane, are studied. Analytical solutions based on the presented uniﬁed model are obtained. Comparisons between results from the presented model and the m(I) model indicate that the explicit partition of frictional and collisional stress components provides insights in dense granular ﬂows. In addition, the new model is used to predict the stress-strain relations in two annular shear tests. The applicability and advantages of the uniﬁed model are discussed. Keywords: constitutive model, granular-ﬂuid ﬂows, solid volume fraction, dynamic 11 viscosity, dense ﬂow regimes 12 Corresponding author. E-mail: wei.wu@boku.ac.at Preprint submitted to Acta Geotechnica July 16, 2020This is a preprint of the paper published in Acta Geotechnica Accepted citation: Guo, X., Peng, C., Wu, W. and Wang, Y., 2020. Unified constitutive model for granular–fluid mixture in quasi-static and dense flow regimes. Acta Geotechnica, pp.1-13. DOI: https://doi.org/10.1007/s11440-020-01044-1 1. Introduction 13 Granular ﬂows are widely presented in industry and nature, e.g., raw material pro- 14 cessing and geophysical ﬂows. With different material properties, external loads and 15 boundary conditions, granular materials can behave like solid, ﬂuid or gas. Debris ﬂows 16 are typical granular ﬂows with viscous interstitial ﬂuid which show solid-like behaviors 17 before failure and ﬂuid-like behaviors after failure. They are normally treated as a ﬂuid 18 continuum with microstructural effect in constitutive modeling for the ﬂuid-like behav- 19 iors [1, 2]. However, a constitutive model that can describe both the solid-like and ﬂuid- 20 like behaviors is required for a complete simulation of the process from quasi-static state 21 to the fast ﬂowing state. In the quasi-static regime, the solid-like granular materials can 22 be sufﬁciently described using plasticity theories in soil mechanics, whereas constitutive 23 equations based on kinetic theory have been proposed for the gas regime [3]. However, 24 most granular ﬂows lie neither in the slowly shearing quasi-static regime, nor in the fully 25 collisional gas regime. In the intermediate ﬂow regime, the granular volume fraction is 26 close to its critical value; thus, the ﬂows can be considered incompressible. This ﬂow 27 regime is usually referred to as dense ﬂow regime in literature. 28 In the dense ﬂow regime, the granular material behaves like a liquid dominated by 29 micromechanics of the granular mass [4]. Two distinct features are observed for dense 30 granular-ﬂuid ﬂows. The ﬁrst is that there exists a threshold yield stress below which 31 the granular media stops ﬂowing; the second is that the granular-ﬂuid ﬂows are shear- 32 rate-dependent. These two features arise from three types of dominant interactions in 33 dense granular-ﬂuid ﬂows: frictional contacts and collisions between solid particles, and 34 viscous drag forces between solid particles and the interstitial ﬂuid. The continuously 35 evolving network of frictional contacts result in the frictional behaviors in the dense 36 granular ﬂows, whereas the viscous rate-dependent behavior has its source in the ﬂuid 37 viscosity and the particle collisions. Which interaction is more dominant depends on 38 several factors, e.g., material property, external load and boundary condition. 39 In a pioneering work on a gravity-free dispersion of solid spheres sheared in Newto- 40 nian liquids, which can be considered as a simple shear test shown in Figure 1, Bagnold 41 2 Figure 1: Schematic of a simple shear test. [5] proposed a constitutive model in which the shear stress has the relation 42 Tv=2:25l3 2mdU dy=k1dU dy(1) in the so-called ’macro-viscous’ regime and 43 Ti=0:042r s(ld)2dU dy2 sina i=k2dU dy2 (2) in the ’grain-inertia’ regime, respectively. In the above expressions, mis the dynamic 44 viscosity of the interstitial ﬂuid; rsanddare the material density and the mean diameter 45 of the grains, respectively; d U=dydenotes the shear strain rate; the tangent of aicorre- 46 sponds to the ratio between the shear and normal stress in the grain-inertia regime; lis 47 the linear concentration deﬁned by 48 l=d s="C0 C1 3 1#1 ; (3) where sis the free distance between two particles, Cis the mean solid volume fraction, 49 andC0is the maximum possible solid volume fraction when l!¥(s=0). In Bagnold’s 50 case, an analytical value for spheres, 0.74, was employed for C0. In general cases, the 51 maximum measured solid volume fractions were found related to the size of the particles 52 and the container dimensions. Usually, C0is replaced by the asymptotic limit of the max- 53 imum measured solid volume fraction C¥as the container dimensions approach inﬁnity 54 3 [6]. 55 A dimensionless value termed ’Bagnold number’ is deﬁned as 56 B=l1 2rsd2(dU=dy) mµTi Tv(4) to characterize the granular ﬂow regimes as macro-viscous (B <40), grain-inertia (B > 57 450), or a transitional state (40 B450) between the two limits. The viscosity of 58 the interstitial ﬂuid and the particle collisions are believed to be dominant factor in the 59 macro-viscous regime and grain-inertia regime, respectively. 60 The normal stress Pis observed proportional to the shear stress Twith a constant 61 ratio, 0.75, in the macro-viscous region which decreases progressively through the tran- 62 sition region, till reaches another constant value at about 0.32 in the grain-inertia region. 63 It is formulated as 64 Tv Pv=tana v=0:75 (5) in the macro-viscous region and 65 Ti Pi=tana i=0:32 (6) in the grain-inertia region, where the tangents of avandaiare the stress ratios in different 66 regimes and related to the material properties. In Bagnold’s experiments, the density of 67 the solid particles is equal to that of the interstitial ﬂuids. The effect of gravity was elimi- 68 nated, and therefore no yield stress was observed in the experiment. This setup highlights 69 the effect of the ﬂuid viscosity and the particle collisions on the ﬂowing behaviors after 70 failure. The models developed based on the experimental observations never take yield 71 criterion into account and cannot be applied to describe the solid-like behaviors before 72 ﬂowing. In addition, the constitutive relations of (1) and (2) are developed for both the 73 speciﬁc ﬂow regimes. The constitutive model for the transition region is unknown, and 74 the two ﬂow regimes are not uniﬁed in a constitutive framework. Bagnold’s model cannot 75 achieve a complete description of the stress state in the process from the quasi-static state 76 to the ﬂowing state. The transition between the two ﬂow regimes cannot be investigated 77 consquently. 78 Another widely accepted constitutive model for dense granular ﬂow regime is the 79 4 m(I)model [7]. It employs an overall phenomenological frictional coefﬁcient related 80 to a dimensionless ’inertial number’. Frictional contact and collision are not explicitly 81 considered in this model. Furthermore, as it is developed for dry granular ﬂow, the drag 82 force imposed on the solid particles by the viscous interstitial ﬂuid is not taken into 83 account. Boyer et al. [8] extended this model to general cases with viscous interstitial 84 ﬂuids. However, the mechanical behaviors inside the yield surface is undeﬁned in both 85 of the models. Thus, the m(I)model cannot be applied to granular materials in the quasi- 86 static regime. 87 This paper aims to develop a framework of constitutive modeling to unify the consti- 88 tutive theories for the solid-like and ﬂuid-like behaviors of a granular-ﬂuid mixture. A 89 simple shearing model is developed based on the framework by employing the Bagnold- 90 type relations for the dense ﬂow regimes and the elastoplastic theory, a Mohr-Coulomb- 91 type relation, for the quasi-static regime. The applicability of the Mohr-Coulomb-type 92 relation to describe the solid-like behaviors is discussed. Based on the effort devoted to 93 examine Bagnold’s pioneer work [9, 10,11,12], the constitutive relations for the macro- 94 viscous and grain-inertia regimes are modiﬁed and uniﬁed in the new model to consider 95 the dynamic behaviors. It was demonstrated that this framework provides a complete 96 description for granular-ﬂuid materials under slow plastic deforming and fast collisional 97 shearing. In this work, dry granular ﬂows in simple conﬁgurations are studied to get an- 98 alytical solutions based on the presented uniﬁed model. The results are compared with 99 the analytical solutions obtained using the m(I)model. In addition, the new model is 100 employed to predict the stress-strain relations of two annular shear tests, the numerical 101 and experimental results are compared. The applicability and advantages of the uniﬁed 102 model are discussed. 103 2. Framework of the constitutive model for granular-ﬂuid ﬂows 104 As pointed out above, the constitutive relations (1) and (2) were developed for two 105 distinct regimes, respectively. The constitutive relation for the transition region is absence 106 so that the continuous description of the ﬂow regimes is impossible. In granular-ﬂuid 107 ﬂows, the intergranular friction, collisions and the ﬂuid viscous force are believed to 108 coexist and transfer momentum simultaneously [13]. Thus, an additive form combining 109 (1) and (2), which reﬂect the viscous ﬂuid ﬂow and the particle collisions respectively, is 110 5 employed for the shear stress of a dense granular-ﬂuid ﬂow. The total rheological shear 111 stress Tris the sum of TvandTi, i.e. 112 Tr=Tv+Ti: (7) Figure 2shows the comparison between experimental data and numerical predictions us- 113 ing Eq. (7). It is found that the direct additive form can correctly capture the overall 114 dynamic behavior and provide an easy yet seamless transition. The material parame- 115 ters introduced in Bagnold’s experiments are used in the stress prediction by combining 116 equations (1) and (2) into (7). 117 Figure 2: Comparison between the shear stress predicted by (7) and experimental data of l=11 in [5]. The constitutive relation for the transition regime is therefore determined. As the 118 shear rate increases, the domination of the mechanical behaviors changes from the ﬂuid 119 viscosity to the particle collisions spontaneously in a dense suspension. Likewise, the 120 total rheological normal stress Prcan be described by 121 Pr=Pv+Pi: (8) Ness & Sun [14] proposed a similar constitutive equation based on the additive form 122 combining viscous effects with dry granular rheology. In their equation, the prolonged 123 frictional contact and the particle collisions are not distinguished. The former is believed 124 to be rate-independent and cannot be sufﬁciently described by rate-dependent rheological 125 6 models. As mentioned before, the yield stress, derived from the prolonged frictional con- 126 tact, was not accounted for in the Bagnold’s model since the effect of the solid particles’ 127 weight was eliminated in the experiments. For general cases of granular-ﬂuid mixture, 128 the contact friction in the quasi-static regime and residual stresses after failure are es- 129 sential parts of a complete model. The constitutive equations for granular ﬂows should 130 satisfy a generalized Mohr-Coulomb type yield criterion as dU=dy!0 [15, 16]. Shibata 131 & Mei [11] assumed the following relation for the quasi-static regime 132 T0=f P0sgndU dy ; (9) where T0andP0are the shear and normal stresses in the quasi-static state; fis an empir- 133 ical coefﬁcient of dynamic friction. The ’sgn’ function is applied to return the sign of the 134 shear strain rate. 135 The particle weight induced stresses are combined with the rheological stresses by 136 following the additive form to yield a constitutive framework as 137 P=P0+Pr=P0+Pv+Pi; (10a) T=T0+Tr=T0+Tv+Ti: (10b) where P0andT0are regarded as the static portion of the framework. TrandPrare the 138 dynamic portion. The framework is akin to the assumption of constitutive modeling for 139 granular materials in [17], in which the stress tensor can be decomposed into a static 140 (rate-independent) part and a dynamical (rate-dependent) part from the perspective of 141 statics. The framework based on the assumption was preliminarily validated in cases of 142 dry granular ﬂow where the viscous terms PvandTvare negligible [18, 19]. The constitu- 143 tive model based on this framework is capable to describe the stress state throughout the 144 shearing process from quasi-static to high-speed shearing stage. 145 The framework can predict steady uniform ﬂows over a range of bed slopes rather 146 than only for one critical slope for dry granular ﬂows on an inclined plane, which is 147 consistent with the experimental observations [16, 20]. 148 7 3. Constitutive model for dense granular-ﬂuid ﬂows 149 In this section, a concrete simple-shearing model for granular-ﬂuid ﬂow is developed 150 based on the framework of (10). A conventional approach is employed to determine the 151 static stresses P0andT0. The speciﬁc expression for the dynamic portion is determined by 152 modifying the Bagnold’s models. Both of the static and dynamic portions are formulated 153 for the solid phase and therefore regarded as effective stress deﬁned in soil mechanics. 154 It is clear that the stresses determined by (1) and (2) will vanish when the solid volume 155 fraction is equal to zero. 156 3.1. Constitutive relation for quasi-static state 157 As indicated by equation (9), the yield stress T0is proportional to the normal stress P0158 with a constant friction coefﬁcient. The empirical friction coefﬁcient, f, is determined as 159 the tangent of the residual friction angle fsince the yield stress is assumed equal to the 160 residual stress in the Mohr-Coulomb type relation. The friction angle can be obtained in 161 simple shear tests [21], and is assumed independent of the shear rate. From (9), the shear 162 stress in quasi-static stage is expressed as 163 T0=P0sgndU dy tanf: (11) Once P0is determined, the shear stress T0will be obtained consequently. The normal 164 stress P0usually stems from the gravity of the solid particles. 165 Figure 3: Schematic of a steady granular-ﬂuid ﬂow on an inclined plane. Take the free surface granular-ﬂuid ﬂow shown in Figure 3as an example, the normal 166 8 stress is the component of gravity perpendicular to the ﬂow plane, 167 P0= (rsrf)Cghcosq ; (12) where rfis the density of the ﬂuid. It is worth to mention that the normal stress deter- 168 mined by (12) is considered as the initial value of P0. The static stress P0will decrease 169 with the increase of shear rate in the dense ﬂow regime since more ﬁerce particle colli- 170 sions lead to less prolonged contact. 171 The simple relation (11) is widely employed as it is easy to be understood and im- 172 plemented in numerical calculations. The expression (12) was proved applicable in the 173 cases without shear softening (or liquefaction in saturated granular material) in the quasi- 174 static state [22]. For the cases with excess pore pressure, more complicated constitutive 175 theories are required for the quasi-static regime. It is stated that granular materials show 176 rate-independent hypoplastic behaviours under monotonous deformations [23]. 177 3.2. Constitutive relations for dense ﬂow regime 178 From (1) and (2), together with (4), two dimensionless quantities are obtained as 179 Tvrsd2 l m2=2:25B (13) and 180 Tirsd2 l m2= (0:042sin ai)B2; (14) which depend only on the Bagnold number B. These relationships cannot describe the 181 experimental observations when the linear concentration, l, is greater than 12 [5, 12]. It 182 implies that the dimensionless quantities are also dependent on the linear concentration. 183 Thus, the dimensionless quantities are extended to the following general expressions as 184 Tvrsd2 l m2=f1(l)B (15) and 185 Tirsd2 l m2=f2(l)B2; (16) where f1(l)andf2(l)are functions of lto be speciﬁed. 186 9 (a) Shear stress in the macro-viscous regime 187 The shear stress in the macro-viscous regime is ﬁrst investigated with the goal to 188 ﬁnd a constitutive model applicable in the whole spectrum of lstudied in Bagnold’s 189 experiments. For the viscous regime, Bagnold derived the expression of total shear stress 190 as 191 Tv= (1+l) 1+1 2f(l) mdU dy; (17) in which Tvis the total shear stress consisting of the contributions of grain and ﬂuid, 192 Tvandtv;f(l)is a function of the linear concentration and reﬂects the amplitude of the 193 shear velocity ﬂuctuation, and f(0) = 0. In order to determine the shear stress of the solid 194 phase, the total shear stress needs to be divided into Tvandtv. Bagnold assumed a simple 195 relation for the ﬂuid contribution in the total shear stress in the grain-inertia regime as 196 [12,24] 197 ti=t0 (1+l); (18) where t0is the shear stress in pure ﬂuid and formulated as mdU=dy. The ﬂuid contribu- 198 tion decreases as the linear concentration increases. This relation is assumed to hold in 199 the macro-viscous regime as well. Thus, 200 tv=m (1+l)dU dy: (19) Based on (17) and (19), the contribution of gain is derived as 201 Tv=Tvtv = (1+l) 1+1 2f(l) 1 1+l mdU dy:(20) Bagnold employed f(l) =lin (17). Thus, the expression of Tvis preliminarily deter- 202 mined as 203 Tv= (1+l) 1+1 2l 1 1+l mdU dy: (21) 10 The dimensionless quantity derived from (21) is 204 205 Tvrsd2 l m2= l3 2+3 2l1 2+1 2l1 2+l3 2 1+l! B; (22) which demonstrates a concrete formula of the proposition of (15). The model (21) is 206 capable to predict the experimental results for granular-ﬂuid mixture with high concen- 207 tration to some extent. For exact prediction, a further modiﬁcation is needed. 208 Figure 4: Schematic of a simple-shearing granular-ﬂuid ﬂow with a stagnant zone. A well-known phenomenon in granular-ﬂuid ﬂows is that the shear stress is relatively 209 insensitive to the solid volume fraction Cwhen the value of Cis below approximately 0.5, 210 but increases rapidly when the solid volume fraction exceeds this critical value [6]. Some 211 researchers [25, 26] used a power series of the volume fraction (or the linear concentra- 212 tion) combined with an exponential term to describe the dynamic viscosity of granular- 213 ﬂuid mixture. These expressions were determined by ﬁtting shear stress versus shear 214 strain rate curves. The mechanism of the phenomenon needs further clariﬁcation. As 215 pointed out by Bagnold [5] and Hanes & Inman [6], the solid volume fraction must be 216 less than a critical value Ccto assure a full shearing to occur. The critical value was found 217 to depend on the packing pattern and lay between 0.53 and 0.65 [11, 24]. For the general 218 case of a gravity ﬂow shown in Figure 4, a stagnant zone with thickness of Hhwill 219 arise at the bottom of the specimen when the mean solid volume fraction Cis greater than 220 Cc. If the total thickness Hrather than that of the ﬂowing zone his used in the calcula- 221 tion of the shear strain rate, an underestimated result will be obtained. In a steady ﬂow 222 11 with C>Cc, the exchange of particles between the ﬂowing and the stagnant zone reaches 223 a balance. The mean solid volume fraction in the upper layer, Cf, is deduced to be un- 224 changed as the critical value Cc. The mean value of the stagnant zone Csshould be greater 225 than Ccto keep static. Otherwise, solid particles will move from the stagnant zone to the 226 ﬂowing zone till Cfreaches the value of Cc. Thus, the experimental observation that the 227 shear stress increases dramatically when the solid volume fraction exceeds the critical 228 value Ccis attributed to both the increase of the dynamic viscosity and the decrease of 229 the thickness of the ﬂowing zone. Therefore, the critical value of solid volume fraction 230 should be taken into account in the simple-shearing model for dense granular-ﬂuid ﬂow. 231 Inspired by the way to take the critical value into account in [11], the following model is 232 proposed for the macro-viscous regime. 233 Tv=K1dU dy; (23) where 234 K1= (1+l) 1+1 2l 1 1+l m 1C Ccn (24) is called effective viscosity; nis a ﬁtting parameter. In Shibata & Mei’s work [11], it was 235 chosen n=1. 236 (b) Shear stress in the grain-inertia regime 237 The constitutive equation for the grain-inertia regime is slightly modiﬁed based on 238 the Equation (16). The experimental results in [5] show that shear stresses in different 239 regimes deviate as a whole with the change of the linear concentration in the logarithm 240 coordinates of rsd2(dU=dy)2versus grain shear stress. It is assumed that the shear stress 241 Tichanges in the same rate with Tvwhen the value of lvaries. Based on this assumption, 242 the original equation for the ’grain-inertia’ regime, (2), is modiﬁed to be 243 Ti=K2dU dy2 ; (25) where 244 K2=0:042R vrs(ld)2sina i (26) 12 and 245 Rv=K1 k1=h (1+l) 1+1 2l 1 1+li 1C Ccn 2:25l3 2(27) is termed as correction factor. 246 3.3. A simple-shearing model for dense granular-ﬂuid mixture 247 Substituting (23) and (25) into (7) yields the expression of shear stress in the dense 248 ﬂow regime 249 Tr=K1dU dy+K2dU dy2 : (28) Based on (5), (6) and (8), the normal stress in the ﬂowing state is expressed as 250 Pr=K1 tanavdU dy+K2 tana idU dy2 : (29) The modiﬁed models are capable of describing the mechanical behaviors of a dense 251 granular-ﬂuid mixture from the ’macro-viscous’ regime to the ’grain-inertia’ regime. As 252 shown in Figure 5, by choosing n=0:2 and Cc=0:65 [6], (28) and (29) can predict 253 Bagnold’s experimental results with satisfactory agreement. 254 Furthermore, based on the framework (10), the dynamic stresses (28) and (29) are 255 combined with the static stresses T0andP0, respectively, to obtain a simple-shearing 256 model as 257 P=P0+K1 tanavdU dy+K2 tanaidU dy2 ; (30a) T=P0sgndU dy tanf+K1dU dy+K2dU dy2 : (30b) This is a simple example of the uniﬁed model for granular-ﬂuid mixture in the quasi-static 258 and dense ﬂow regimes. The rate-independent constitutive component for the uniﬁed 259 model can be chosen according to the requirement in comprehensively modeling the 260 behaviors inside the yield surface. 261 13 (a) (b) Figure 5: Comparison between the stresses predicted by (28) and (29) with n=0:2 and Cc=0:65 and the data from (a) the ﬁgure 3 and (b) the ﬁgure 4 in [5]. The experimental data are indicated by various symbols. The dashed lines (upper panel) denote the shear stresses while the solid lines (lower panel) are the normal stresses. 14 4. Analytical solutions for steady granular ﬂows 262 Three dry granular ﬂows under simple conﬁgurations are studied. Analytical solu- 263 tions from the simple-shearing model are compared to the solutions from the m(I)model. 264 4.1. Plane shear ﬂow 265 Consider shearing of a granular material made of particles with diameter dand den- 266 sityrs, as shown in Figure 1. The conﬁning pressure Pis imposed on the top plate. The 267 granular material is sheared at a prescribed rate ˙g=U=h, where Uis the velocity of the 268 top plate and his the depth of the granular body. Gravity is not considered. To maintain 269 a steady plane shear, the shear force on the top plate depends on the imposed conﬁning 270 pressure and shear rate. 271 In the m(I)model, the inertial number can be calculated through 272 I=d˙gp P=r s: (31) Consequently, we can obtain the shear stress Tby [7] 273 T=m(I)P;with m(I) =m1+m2m1 I0=I+1; (32) where m1,m2andI0are three parameters of the m(I)model. 274 With the presented uniﬁed model, considering the steady shearing of dry granular 275 materials where the viscous terms are negligible, we have the following two expressions 276 T=sf yysinf+K2˙g2; (33) and 277 P=sf yy+K2 tanai˙g2: (34) The frictional shear stress is computed as sf yysinf, instead of sf yytanf, because sim- 278 ulations of plane shear ﬂows indicate that the difference between the vertical and the 279 horizontal normal stress is very small, i.e., sf xxsf zz[28,29]. The frictional shear stress 280 15 is then equal to sf yysinf. Based on Eqs. (33) and (34), the total shear force is obtained as 281 T=Psinf+K2˙g2(1sinf tanai): (35) Comparing Eqs. (32) and (35), it can be observed that the shear stresses obtained by 282 using the m(I)model and the uniﬁed model have a similar structure. Both results consist 283 of two parts, one part is constant with given normal stress P, and the other one depends on 284 the shear rate. In the dense ﬂow regime, the shear stress obtained using the m(I)model 285 ranges from m1Ptom2P. On the other hand, in the dense ﬂow regime, the dispersive 286 pressure Pi=K2˙g2=tanaiis from 0 to P; therefore, the range of shear stress obtained 287 using the uniﬁed model is from PsinftoPtanai. It is observed that the constants sin f 288 and tan aiin the uniﬁed model are equivalent to the constants m1andm2in the m(I) 289 model. 290 Both of the models enable precise classiﬁcations of the three ﬂow regimes of dry 291 granular materials: quasi-static, dense ﬂow and gaseous. Small values of ˙g(˙g!0) 292 correspond to the quasi-static regime in the sense that deformation is very slow. In the 293 m(I)model, the threshold of Icorresponds to the rapid and dilute ﬂow regime is I1, 294 which is obtained from experiments [29]. Contrarily, in the uniﬁed model, the gaseous 295 regime can be theoretically determined: it corresponds to the state when the network 296 of frictional contacts disappears, that is, when sf yy=0. Based on the aforementioned 297 analysis, the shear rate corresponds to the dense ﬂow regime in the plane shear is 298 for the m(I)model : 0 <˙g.p P=q rsd2; for the uniﬁed model : 0 <˙g<p P=pK2:(36) Again, the ranges of shear rate corresponding to the dense ﬂow regime obtained using 299 them(I)model and the uniﬁed model have similar form. Note that the constant K2is 300 related to particle diameter, density and solid volume fraction. It can be observed that in 301 both of the models, the upper limits of the shear rate depend on the conﬁning pressure P 302 and physical properties of the granular material. 303 Although these two models bear some similarities, there are fundamental differences. 304 Them(I)model employs an inertial-number-dependent phenomenological frictional co- 305 16 efﬁcient, which changes with shear rate and conﬁning pressure. However, the model does 306 not give any clue why the phenomenological frictional coefﬁcient m(I)changes with dif- 307 ferent shear rate and conﬁning pressure. On the other hand, the uniﬁed model explicitly 308 considers the origin of different stress components. In the uniﬁed model, the change 309 of phenomenological frictional coefﬁcient can ﬁnd its explanation as follows. The ratio 310 between shear stress and normal stress in the network of frictional contacts is related to 311 the internal friction angle fof the material, whereas the ratio of the collision-induced 312 shear/normal stresses is related to another angle ai. Experiments show that aiis gener- 313 ally larger than f. Therefore, when the ﬂow state moves from the quasi-static regime to 314 the gaseous regime, the collision-induced stresses becomes more and more signiﬁcant. 315 In this process, the phenomenological frictional coefﬁcient increases from sin fto tan ai. 316 The evolution of the phenomenological frictional coefﬁcient in the two models is written, 317 respectively, as 318 in the m(I)model : m=m1+˙gd(m1m2) I0p P=r s+˙gd; in the uniﬁed model : m=sinf+˙g2K2(tanaisinf) Ptanai:(37) Comparing Eq. (37) with Eq. (32), the evolutions of the phenomenological frictional 319 coefﬁcient in the two models have a similar structure, but different dependences on shear 320 rate and conﬁning pressure. 321 4.2. Vertical chute ﬂow 322 Consider a granular material ﬂowing between two vertical plates with a distance of 323 2W, as shown in Figure 6. The conﬁning pressure on the two vertical plates is P, which 324 implies that sxx=Pwithin the ﬂow in the steady state. 325 According to the uniﬁed model, the possible range of shear stress txyin the chute ﬂow 326 is from PsinftoPtana. Because of symmetry, the shear stress at the center of the ﬂow 327 is zero. Based on the force balance along the xaxis, the shear stress can be written as 328 txy=rgx; (38) 17 Figure 6: Schematic of vertical chute ﬂow. where ris the bulk density of the granular material. The steady ﬂow condition requires 329 that the shear stress increases linearly from the ﬂow center to the vertical plates. Con- 330 sequently, there exists an unsheared zone in which txy<Psinf. Namely, when x<xp 331 (xp=Psinf=rg), the granular material is unyielded and moves as a plug layer. If xp>W, 332 the whole granular body is unyielded and no granular ﬂow is possible. Furthermore, the 333 maximum shear stress required for the force balance does not exceed the maximum shear 334 stress that can be provided, i.e. txy<Ptanai. Therefore, according to the uniﬁed model, 335 a steady granular ﬂow is possible in the vertical chute only if 336 rgW=tanai<P<rgW=sinf: (39) Now we proceed to derive the velocity proﬁle of the steady ﬂow. The stress compo- 337 nents for the present vertical chute ﬂow, as speciﬁed in (33) and (34) for the plane shear 338 ﬂow, take the forms 339 P=sf xx+K2 tanaidUy dx2 ; txy=sf xxsinf+K2dUy dx2 :(40) Substituting Eq. (40) into Eq. (38), we have 340 dUy dx2 =(rgxPsinf)tanai K2(tanaisinf): (41) Solving the above equation using the boundary condition Uy=0 atx=W, we can obtain 341 18 the velocity proﬁle 342 Uy(x) =Up y+2ptanai(rgxPsinf)3=2 3p K2(tanaisinf)rgifxpxW; (42) where Up yis the velocity of the unsheared plug layer 343 Up y=2ptanai(rgWPsinf)3=2 3p K2(tanaisinf)rgifx<xp: (43) Analytical solutions can also be obtained using the m(I)model [30]. Similarly, an 344 unsheared plug layer exists, and there is a conﬁning pressure range, only with conﬁning 345 pressure in this range a steady granular ﬂow is possible. If the parameters in the m(I) 346 model are taken as m1=sinfandm2=tanai, the width of the plug layer and the range 347 of the conﬁning pressure are identical to those obtained from the uniﬁed model. However, 348 the velocity proﬁles obtained using these two models are not identical. 349 To quantitatively compare the solutions obtained using these two models, the velocity 350 proﬁles are calculated. The following material properties are used: r=1470 kg/m3, 351 rs=2450 kg/m3andd=0:002 m. The parameters for the m(I)model are chosen as 352 m1=0:38, m2=0:64 and I0=0:279. For the uniﬁed model, we take sin f=m1, tana= 353 m2andK2=0:2. The width of the chute is 0.2 m, the conﬁning pressure Pis 3000 354 Pa. Owing to the symmetric geometry, only the right half of the granular material is 355 considered. The calculated results are shown in Figure 7. Using the given parameters, 356 the two models predict identical unsheared velocity. However, the velocity proﬁles in the 357 shear zone are somewhat different. 358 4.3. Flow on an inclined plane 359 The problem of a granular mass on an incline is considered. As shown in Figure 3, 360 the granular mass has a free surface and is subjected to gravity. The inclination angle is 361 q. The height of the granular layer is h. 362 For this steady granular ﬂow, the stress distribution should be hydrostatic 363 syy=rg(hy)cosq;sxy=rg(hy)sinq: (44) 19 Figure 7: Velocity proﬁles obtained using the uniﬁed model and the m(I)model in the vertical chute ﬂow. For an inﬁnite inclined plane, we assume that the frictional stresses sf yy=sf xxin the 364 ﬂowing state; thus, the shear/normal stress ratio of the frictional contacts is tf xy=sf yy= 365 sinf. 366 The following force balance must be satisﬁed in the steady state 367 syy=rg(hy)cosq=sf yy+K2 tanaidU dy2 ; txy=rg(hy)sinq=sf yysinf+K2dU dy2 :(45) Here, once again, the stress consists of the contributions from frictional contacts and 368 collisions, described by the Mohr-Coulomb model and the Bagnold-type rheology, re- 369 spectively. Based on the force balance equations, the following differential equation is 370 derived 371 (1sinf=tanai)K2dU dy2 rg(sin qcosqsinf)(hy) = 0: (46) Solving the above equation, the velocity proﬁle is obtained as 372 U(y) =2 3s rg(sin qcosqsinf) (1sinf=tanai)K2(h3=2(hy)3=2): (47) 20 The velocity proﬁle, evolving as (hy)3=2, corresponds to the Bagnold proﬁle. Once the 373 velocity is obtained, the stress components from frictional contacts and collisions can be 374 determined. From the velocity proﬁle it is noted that velocity and shear rate increase with 375 the inclination q. With increasing q, the normal and shear stress induced by collisions 376 become more signiﬁcant. This increase in the dispersive normal stress results in a trend 377 of dilation, tending to reduce the stress components caused by frictional contacts. 378 If the granular material is modeled using the m(I)model, dense granular ﬂow can only 379 happen if the inclination qis between tan1m1and tan1m2. With the uniﬁed model, a 380 similar inclination range can be derived as follows 381 tanq=txy syy=sinf+sc zz szz(tanaisinf); (48) where sc zzis the dispersive normal stress. The steady dense ﬂow condition requires 382 0<sc zz<szz, because if sc zz>szzthe network of frictional contacts disappears com- 383 pletely, which means that the granular material enters the gaseous regime. Consequently, 384 according to the uniﬁed model, the inclination range for dense granular ﬂow is 385 tan1(sinf)<q<ai; (49) which is exactly the inclination range derived from the m(I)model. It is easily seen that 386 the obtained velocity proﬁle is valid only under the condition (49). 387 The velocity proﬁles of the granular ﬂow with two different inclinations are cal- 388 culated. The material parameters used in the calculation are the same as those used 389 in the former vertical chute ﬂow. For the uniﬁed model, the parameters are taken as 390 sinf=m1=0:38 and tan ai=m2=0:64. Figure 8shows the results obtained using 391 K2=0:3 for q=23, and K2=0:1 for q=25. Using the given parameters, the two 392 models provide similar predictions. 393 5. Performance of the proposed model 394 In this section, the performance of the uniﬁed simple shear model is demonstrated 395 by modeling two tests, respectively, for a dry granular ﬂow and a granular-ﬂuid mixtures 396 ﬂow. 397 21 Figure 8: Velocity proﬁles obtained using the uniﬁed model and the m(I)model in the ﬂow on inclined plate. 5.1. Case 1: Dry granular materials 398 The experimental results of dry granular materials sheared in an annular shear cell 399 were reported by Savage & Sayed [21]. In the experiments, the shear velocity was 400 adjusted to keep the height of the granular mass constant and thus keep the volume 401 unchanged. This steady state is equivalent to an undrained simple shear of saturated 402 granular materials. The experimental data for 1.0 mm spherical polystyrene beads are 403 predicted by the new model. The shear stresses were measured on the top surface of the 404 sample. As stated in [21], the loads applied by the upper disk range from 100 to 1500 405 N=m2. However, the speciﬁc load for each single sample was not reported. By checking 406 the measured normal stress for 1.0 mm beads, we assume that the normal stress in the 407 quasi-static regime P0has a value around 100 N=m2which is also regarded as the resid- 408 ual normal stress after yielding. As mentioned before, the theoretical maximum value 409 of solid volume fraction, C0, will be replaced by the asymptotic limit of the maximum 410 measured solid volume fraction, C¥, in the calculation for a speciﬁc experiment. Since 411 the exact value of C¥was not reported in Savage & Sayed’s work, we use 0.65 which is 412 a typical value for monosized spheres [6, 24]. The critical volume fraction Ccis approxi- 413 mately 0.62 [11]. All the parameters for the prediction are listed in Table 1and Table 2. 414 415 As shown in Figure 9, the predicted curves can ﬁt the experimental data at high 416 22 Table 1: Parameters for the shear tests of 1.0 mm beads in [21] d C ¥ Cc rs rf m f q [mm] [-] [-] [kg=m3] [kg=m3] [Pas] [] [] 1.0 0.65 0.62 1095 1.29 1:83 10523 0 Table 2: Stress ratios measured in the shear tests of 1.0 mm beads in [21] C tanav tanai [-] [-] [-] 0.461 0.50 0.51 0.483 0.40 0.51 0.504 0.30 0.51 0.524 0.30 0.51 shear strain rate very well. For the dense specimen, the non-quadratic dependence of 417 the stresses on the shear rate in the slow shear stage is also captured by the new model. 418 However, the predicted results for the specimens with C<0:524 show a slight over- 419 estimation, especially in the slow shear stage. The assumed value of P0, 100 N=m2, 420 is believed to be greater than the real value of the tests with C<0:524. However, as 421 mentioned before, P0=100 N=m2is the minimum value which can be imposed on the 422 specimens by the upper disk. A plausible explanation is that shear softening occurs in the 423 relatively loose specimens under undrained condition. In this case, the normal stress P0424 which corresponds to the effective stress in soil mechanics would decrease to a residual 425 value to eliminate the tendency of volume compression. A looser specimen will present 426 lower residual strength. The Mohr-Coulomb type relations (11) and (12), used as the 427 static portion of (30), cannot capture the decrease of P0in the quasi-static stage. How- 428 ever, the dynamic portion is validated in the describing of the stress-strain rate relations 429 for dry granular ﬂows with different solid volume fraction. 430 23 (a) (b) Figure 9: Elementary modeling results for dry granular ﬂow with different grain linear concentration: (a) normal stress and (b) shear stress. The experimental data are indicated by various symbols. The dashed lines (upper panel) denote the normal stresses while the solid lines (lower panel) are the shear stresses, predicted by the postulated model 24 5.2. Case 2: Granular-water mixture 431 Another elementary modeling by the uniﬁed simple-shearing model is based on the 432 laboratory tests from Hanes & Inman [6] with spherical particles sheared in water. The 433 data for particles with diameter of 1.85 mm are selected since they are claimed to be 434 of good quality. The maximum measured volume fraction for 1.85 mm particles was 435 reported to be 0.55 [6]. The asymptotic limit C¥is presumed to be approximately 0.61. 436 The critical volume fraction is assumed to be 0.52 as a partially shearing with C=0:53 437 was reported. The load from the upper disk, P0, is about 200 N=m2. The angle avis equal 438 to the dynamic angle of repose, 28. All the parameters are listed in Table 3and Table 4. 439 Table 3: Parameters for the shear tests of 1.85 mm beads in [6] d C ¥ Cc rs rf m f q [mm] [-] [-] [kg=m3] [kg=m3] [Pas] [] [] 1.85 0.61 0.52 2780 1000 1:0 1030.59 0 Table 4: Stress ratios measured in the shear tests of 1.85 mm beads in [6] C tanav tanai [-] [-] [-] 0.49 0.53 0.59 0.51 0.53 0.59 440 The predicted results are shown in Figure 10. The predicted curves of C=0:51 441 show a good agreement with the experimental data. The stresses of the specimen with 442 C=0:49 are slightly overestimated. Similar to the case 1, the overestimation of P0and 443 T0may be responsible for this discrepancy. Nevertheless, the dynamic portion of the 444 presented model is capable of describing the stress versus shear rate relation in the dense 445 ﬂow regime for a granular-water mixture. 446 25 (a) (b) Figure 10: Elementary modeling results for granular-water ﬂow with different solid volume fraction: (a) normal stress and (b) shear stress. The experimental data are indicated by various symbols. The solid lines denote the normal stresses in the upper panel and the shear stresses in the lower panel, predicted by the postulated model. 26 6. Conclusions 447 A framework for the constitutive modeling of granular-ﬂuid mixture is developed by 448 combining a rate-independent static part with a rate-dependent rheological part in an ad- 449 ditive form. In this framework, the effect of the prolonged contact, the viscosity of the 450 interstitial ﬂuid and the particle collisions are considered explicitly. The Mohr-Coulomb- 451 type relation for the quasi-static regime and the Bagnold-type rheological model for the 452 dense ﬂow regime are uniﬁed in the framework to obtain a simple-shearing model which 453 can describe the mechanical behaviors throughout the shear process from yielding to 454 high-speed shearing. In the modiﬁed Bagnold-type model, the conversion between the 455 so-called macro-viscous and grain-inertia regimes can be captured continously. The criti- 456 cal solid volume fraction to assure the occurrence of a full shearing is taken into account. 457 The existence of the stagnant zone is identiﬁed as the main reason for the dramatic in- 458 crease of the dynamic viscosity as the solid volume fraction exceeds the critical value. 459 Analytical solutions for dry granular ﬂows in three conﬁgurations, i.e., plane shear, 460 vertical chute ﬂow and ﬂow on inclined plate, are derived based on the presented model. 461 The solutions are compared with those obtained using the m(I)model. The following ob- 462 servations can be made. (1) Both models show that in order to switch from the quasi-static 463 to the gaseous regime, one can either increase the shear rate or decrease the pressure. 464 (2) The overall phenomenological frictional coefﬁcient is dependent on ﬂow kinematics, 465 material properties and boundary conditions. In both models there are lower and upper 466 limits for this phenomenological frictional coefﬁcient. The m(I)model shows that the 467 friction angle increases as the granular material approaches the gaseous regime, but does 468 not give a physical background. The presented uniﬁed model shows that the increase in 469 the phenomenological frictional coefﬁcient is because of particle collisions. 470 The developed model is used to predict the stress-strain rate relations of two granular- 471 ﬂuid ﬂows. The predicted results for specimens with different solid volume fraction show 472 a good agreement with the experimental data in the high shear rate stage. However, when 473 the static (or residual) normal stress in the beginning of the ﬂowing is determined as the 474 component of gravity or the external load on the solid phase, the total stresses in the low 475 shear rate stage is overestimated by the new model for the relatively loose specimens. A 476 possible explanation is that shear softening occurs in the relatively loose specimens in 477 27 the quasi-static regime and the Mohr-Coulomb model cannot capture this behavior. 478 Nevertheless, the uniﬁed model has some favourable distinct features. First, the quasi- 479 static solid-like behaviors of granular materials are fully deﬁned in the rate-independent 480 part of the model. The modeling of the quasi-static regime, such as critical state, plas- 481 tic ﬂow and nonlinear elasticity, can be further improved by changing the static part in 482 the uniﬁed model to more sophisticated models. Thus, the model is quite ﬂexible for a 483 modiﬁcation. Second, granular materials in the quasi-static and dense ﬂow regime can be 484 described using the model uniformly, which simpliﬁes numerical implementation. Third, 485 the model reveals that in the dense ﬂow regime, especially when granular materials are 486 quite close to the gaseous regime, particle collisions contribute to a signiﬁcant part of the 487 stress. 488 Acknowledgements 489 The authors wish to thank the ﬁnancial support from the National Natural Science 490 Foundation of China (grant number 51809230 and 51709230), the European Commis- 491 sion 7th Framework Programme project MUMOLADE (289911), and the H2020 Marie 492 Skłodowska-Curie Actions RISE 2017 HERCULES (778360) and FRAMED (734485). 493 Conﬂict of interest: 494 The authors declare that they have no conﬂict of interest. 495