Transient Evolution of Rheological Properties of Dense Granular Inertial Flow Under Plane Shear

Exploration of the transient evolution of the rheological properties of dense inertial flow can reveal the equilibrium mechanism of granular materials maintaining their own stability under shear. Here, discrete element method simulations are performed to study the transient flow characteristics of a dense granular system under plane shear in the inertial regime. We quantitatively analyze the changes in the system’s flow state, interfacial friction coefficient, effective friction coefficient, microstructure anisotropy, and internal shear strength. Simulation results show that the evolution of the horizontal flow experiences three typical stages, namely transmission, adjustment, and stabilization. Moreover, the shear dilatancy caused by the vertical movement of particles, gradually weakens the spatial geometric constraint and reduces the system’s tangential load-bearing capacity, thereby decreasing the interfacial friction coefficient μ′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mu }^{^{\prime}}$$\end{document}, which represents the boundary driving strength. On the other hand, the shear flow induces variations in the anisotropies of both contact orientation and contact forces, thus improving the system’s internal shear strength μ∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mu }^{*}$$\end{document} and increasing the effective friction coefficient μe\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mu }_\text{e}$$\end{document}. By comparison, μ′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mu }^{^{\prime}}$$\end{document} is greater than μe\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mu }_\text{e}$$\end{document} until approximately equal in the steady flow state. Therefore, during the evolution of the flow state, the boundary driving strength is reduced while the system’s shear resistance is enhanced, and eventually a balance between them is achieved. Distinguished from the micromechanical behaviors, the internal shear strength always mainly originates from the anisotropies in contact orientation and in normal contact force. Moreover, the contribution of the anisotropy in contact orientation becomes more predominant with increasing shear velocity.


Introduction
Under shear, granular materials exhibit fascinating properties such as sustaining stress as a solid or flowing as a fluid [1][2][3]. Granular flow has motivated many studies over the past decades, since the rheological properties of dense granular materials play an important role in various natural phenomena and industrial applications such as geographical faults, powder metallurgy, three-body friction, and particle flow lubrication [4][5][6][7][8][9][10]. Generally, three different granular flow regimes are categorized according to the flow rate, including the quasi-static regime, inertial regime, and collisional regime [11][12][13][14]. In the inertial regime, when a dense granular system is subjected to continuous shear, a local flow near the boundary is first generated after the initial static equilibrium is broken, and gradually a global flow is achieved and ultimately maintained at steady state [15,16]. Moreover, the evolution of the flow state is accompanied by dilatancy, which represents the volume expansion of the granular system [17][18][19][20]. Two representative quantities, i.e., the macroscopic friction coefficient and the solid fraction , are most commonly employed to characterize the rheological properties of dense granular flow [21][22][23][24][25]. Depending on the calculation method, the macroscopic friction coefficient can be referred to three different meanings, namely, the interfacial friction coefficient ′ , the effective friction coefficient e , and the internal shear strength * . The interfacial friction coefficient ′ can be understood as the driving strength provided by the boundary for maintaining the particle flow, which is defined as the ratio of the shear stress to the normal stress on the shear boundary [23][24][25][26][27][28]. The effective friction coefficient e and the internal shear strength * represent the shear resistance to the flow of the system. When the granular system is assumed to be a continuous medium for stress analysis [29], the effective friction coefficient e is expressed as the ratio of the deviatoric stress to the average stress in the system [30,31]. The internal shear strength * is derived from the anisotropies of the microstructure formed by contacts [21,32]. In the steady state, the boundary driving strength and the system's shear resistance are in equilibrium, and these three coefficients ′ , e and * are approximately equal and can be unified as the macroscopic friction coefficient .
The interfacial friction characteristics could be influenced by several geometric factors, including the boundary surface structure, the thickness of particle layers, and the shape and constrained extent of particles [33][34][35][36][37][38]. Most experimental studies on the interfacial friction characteristics address local shear deformation, where stick-slip or particle flow occurs only near the boundary due to either the small shear velocity or the short shear strain. However, when the boundary cumulative shear strain is large enough in the inertial regime, the system reaches a steady flow state eventually. In this case, the interfacial friction coefficient ′ tends to be a constant value independent on the initial filling state of the dense granular system. Moreover, the interfacial friction coefficient ′ has a significant dependence on the flow rate, which is determined by the shear velocity and the applied normal stress [39][40][41]. For rigid dry granular materials, the flow rate can be measured by the dimensionless inertial number I, which represents the ratio of the inertial force to the confining pressure [14,42,43]. Numerous experimental and simulation results have confirmed that in the inertial regime, the interfacial friction coefficient ′ monotonically increases as a function of the inertial number I [23][24][25][43][44][45]. It is a familiar feature of granular rheology that the interfacial friction coefficient ′ is closely related to the system's flow rate during the shear process.
On the other hand, the interparticle contact and contact forces exhibit anisotropy, which are the microscopic origins of the system's shear resistance [12,21,[46][47][48]. At the microscopic scale, particles transfer forces through interparticle contacts. The variation in the flow state changes the contact status and the force transmission paths between particles, thereby altering the topology and mechanical properties of the microscopic contact force network [47][48][49][50]. As a result, the internal shear strength significantly changes during inertial flow due to the temporal and spatial variations in the microstructure of the contact force network between particles. Much work has been carried out on analyzing the internal shear strength under static filling or steady inertial flow situations. For the stationary granular system, the internal shear strength is independent of the size polydispersity of the filled particles [51][52][53]. In addition, the influences of the different particle shapes on the internal shear strength were also investigated [12,33,34,38,[54][55][56][57][58]. For the steady inertial flow state, Azéma et al. [12] compared the rheological properties of a frictionless granular system between pentagonal particles and circular discs. Binaree et al. [34] further explored the integrated effects of interparticle friction and particle shape on the internal shear strength. In addition, the rheological properties and microstructure characteristics of some cohesive granular systems were investigated [59][60][61], where the effect of viscosity on the internal shear strength was discussed.
By virtue of the abovementioned studies, the steady rheological properties of dense granular materials have been well studied under equilibrium conditions. Under some circumstances, it is necessary to study the transient behaviors, for example, whether a steady flow can be obtained for a particular granular system, and how long it takes for a static system to achieve the steady flow state. To answer these questions, we should understand how the equilibrium is gradually established, and figure out which factors could affect the evolution process. To this end, it is essential to perform a dedicated investigation into the transient evolution of the rheological properties. Here we choose one of the most commonly studied granular flow models-the plane shear model, as our research target. The composing structure and motion manner of this model are simple. Additionally, under the steady flow state, the flow velocity linearly increases along the height, which can be used as a direct criterion for judging whether the system reaches steady state [62]. The particle flow of this model is still a hot topic very recently [39,59,63]. Discrete element method (DEM) serves as a powerful tool for studying the mechanical properties of granular materials because it can extract the microscopic information of individual particles and also analyze the macroscopic behaviors of the whole granular system [64,65]. Therefore, this paper employs DEM simulations to analyze the evolution flow process of a plane shear system by carefully characterizing the evolution principles of the flow state, the mechanical behaviors, and the microstructure characteristics at different moments. Moreover, we clarify the influences of the flow state on the interfacial friction coefficient, the effective friction coefficient and the internal shear strength and finally reveal the equilibrium mechanism between the boundary driving strength and the system's shear resistance. In addition, the influences of the shear velocity are explored.

DEM Model of Plane Shear
In the two-dimensional DEM model, the granular material is described as a collection of discrete discs. The contact force dependent on the force-displacement relationship between particles is defined by the soft contact method. The motion information of each particle is derived from its own unbalanced force using Newton's second law. As shown in Fig. 1a, the dense granular system consists of N p = 2 × 10 4 cohesionless rigid disc particles with an average diameter of d = 1 mm . The solid fraction of the system is 0.84. Initially, the system's width along x-axis is L = 200d , and the height along y-axis is H = 96d . Compared with similar studies, the number of particles included in this simulation model is large enough to ensure the statistical confidence of the analysis results. The particle size distribution polydispersity is set as 20% to avoid crystallization. The top and bottom plates are rough boundaries formed by the bonding monodisperse particles with diameter d , while the sidewalls are periodic boundaries. During shear, the bottom plate is fixed, and the top plate slides tangentially at a constant velocity V s and applies a constant normal stress P to the granular system.
For cohesionless rigid particles, a linear contact model (spring-dashpot model) is used to define the contact forces between particles, as shown in Fig. 1b. The linear model consists of two parts, namely, a stiffness model and a slip model. The stiffness model provides a linear elastic relationship between the contact force and the relative displacement.
The normal contact force f n c is composed of the normal linear force f n cl and the normal damping force f n cd , written as f n c = f n cl + f n cd = k n g s + d nġs , where k n is the normal contact stiffness, d n is the normal damping coefficient, g s is the surface gap between particles, and ġ s is the relative normal translation velocity. For dry granular materials, the normal contact force between particles is always compressive, so the equation f n cl + f n cd ≤ 0 is satisfied. The tangential contact force f t c is incrementally updated using f t where f ot c is the tangential contact force before updating, k t is the tangential stiffness, Δ s is the relative tangential displacement increment, d t is the tangential damping coefficient, and Δ̇s is the relative tangential translation velocity. The relative motion between particles follows the Coulomb relationship | | f t c | | ≤ c f n c , where c is the sliding friction coefficient between particles, dictating that relative sliding between particles is allowed only when the equal sign is satisfied.

Parameters Setup of DEM Model
A series of dimensionless parameters related to the particle material properties and the shear conditions are set according to references [38,43], which mainly include the contact stiffness coefficient k , the stiffness ratio s , and the restitution coefficient e . The dimensionless contact stiffness coefficient k = k n ∕P measures the rigidity of particles and limits the maximum deformation of particles under external load, where P is the applied normal stress, and k n is the normal contact stiffness. In detail, a large value of k represents rigid particles, and a small value represents soft particles. Relevant studies have shown that once k exceeds 10 4 , it will no longer influence the results because rigid particles are ensured [23,66]. Therefore, k is set to 10 4 herein. The stiffness ratio s = k t ∕k n takes an empirical value of 0.5. The normal and tangential damping coefficients d n and d t are related only to the restitution coefficient e once the stiffness coefficient k n and the stiffness ratio s are determined, that is, d n = 2 n √ mk n and d t = 2 s √ msk n , where s and n are normal and shear critical-damping ratios, respectively, and s = n = − ln(e) 2 + (ln(e)) 2 . The restitution coefficient e determines the threshold value for granular system transformation from dense flow to dilution collision [14]. To broaden the range of inertial flow, e is set as a fixed value of 0.1. The particle-particle and particle-wall friction coefficients are both set to 0.4 and notated as c . A combined set of e = 0.1 and c = 0.4 is frequently used in simulations because it results in a strong dissipative granular system that well matches the real situation. Table 1 lists the detailed parameters of the granular system.
The inertial number, defined as I =̇d∕ √ P∕ , is one of the most important quantities to express constitutive laws [25], where ̇ is the shear rate,d is the mean particle diameter, P is the normal stress, and is the particle density. For a frictional granular material, a small value of I ( ≤ 10 -2 ) corresponds to the quasi-static regime, while a large value of I ( ≥ 10 -1 ) for the collisional regime [22,24]. In this study, as P is a fixed value of 50 Pa, the shear velocity of the top plate is selected within the range of 0.5-2.5 m/s, fully ensuring an inertial flow with I in the range of 0.033 to 0.133. The DEM simulation adopts an explicit numerical scheme in the calculation process, where Newton's second law needs to be repeatedly called to update the particle's position and interparticle contact information. During the simulation process, the instantaneous micromechanical information is obtained by an interval of every 1000 steps, including the particle's position and velocity, contact position and contact forces, to monitor the evolution of the flow state and multiscale mechanical behaviors.

Microscopic Quantities
For a two-dimensional granular system, the solid fraction is defined as = A p ∕A , where A p is the area occupied by the particles and A is the total area of the system. The coordination number Z is the average number of contacts and is calculated using Z = N c ∕N p , where N c is the total number of contacts, and N p is the total number of particles in the system. The existing contacts can be grouped into strong contact and weak contact using the following rule: if the force magnitude f c of a specific contact is greater than the average magnitude of all the contact forces ⟨f c ⟩ , namely, f c ≥ ⟨f c ⟩ , this contact is defined as a strong one; otherwise, if f c < ⟨f c ⟩ , it is regarded as a weak one. The number ratio of strong contacts is defined as R s c = N s c ∕N c , where N s c is the number of strong contacts and N c is the total number of contacts. Additionally, the bearing ratio of strong contacts is calculated using the sum of the strong contact forces divided by the sum of all the contact forces, given as where c ∈ s represents the collection of the strong contacts. Likewise, for the remaining weak contacts, the following equations satisfy R w

Macroscopic Quantities
The strain rate tensor is expressed as the symmetric part of the velocity gradient tensor [67][68][69]: where V and V are the velocity fields. In our model, since we assume that the flow characteristics of particles in the sliding direction x are homogeneous, V x ∕ x and V y ∕ x are always 0. Moreover, as the value of V y is very slight, V x ∕ y composes the main variation in the strain rate tensor.
To measure the transient evolution of the global flow in the granular system, a dimensionless equivalent shear strain rate ̇ε e xx is calculated by averaging the horizontal velocity of particles in the system with area A, given as follows: where N p is the total number of particles, V s is the constant sliding velocity of the top plate, and V i x is the horizontal velocity of particle i. In the evolution process, ̇ε e xx <1.0, and when the system reaches steady state, ̇ε e xx ≈ 1.0. The stress tensor is the sum of static contact stress and kinetic contact stress [68,70,71], given as follows: where A is the area of the model, r c is the contact branch vector, f c is the contact force, m i is the mass of particle i , and V i is the velocity fluctuation of particle i relative to the average velocity of the corresponding layer.

Calculation Methods of Friction Coefficients
For the plane shear system, the interfacial friction coefficient ′ is the ratio of the total tangential stress t to the total normal stress n at all contact points between the top plate and the top layer particles, where t c and n c are the tangential stress and the normal stress on a single contact point, respectively, and c ∈ B represents the collection of contacts between particles and the top plate. The effective friction coefficient e is defined as the ratio of internal shear stress to normal stress, and e = xy ∕ yy , where xy and xy are stress tensor components defined in Eq. (3).

Stress-Force-Fabric Relationship
The stress-force-fabric relationship provides an effective approach for establishing the quantitative relationship between the macroscopic stress tensor and the microstructure anisotropy parameters [72]. Based on this relationship, we conduct a statistical analysis on the distribution functions of contact orientation, normal contact force, and tangential contact force in polar coordinates, denoted as P(θ), f n ( ) , and f t ( ) , respectively, where is the angle of unit vector n along the contact direction. These distribution functions provide abundant micromechanical information for accurately evaluating the anisotropy extent of granular systems. Commonly, low-order Fourier approximations are employed to describe them: where f n is the magnitude of average normal contact force over the system, a c represents the anisotropy in contact orientations, a n and a t describe the directional variation of contact forces, c is the direction of anisotropy in contact orientation, and n and t are the preferred directions of contact forces. When the three main principal directions are assumed to be collinear and codirectional, the internal shear strength * can be approximately given as follows: This equation shows the microscopic components of the internal shear strength of the granular system under plane shear.

Transient Evolution of the Flow State
In this part, we employ the simulation results of P = 50 Pa and V s = 0.5 m/s as a demo to quantitatively analyze the evolution of flow state. The simulation time is 40 s, which is long enough for the system to reach steady state. The flow state is characterized by the flow features of both horizontal and vertical directions. In detail, the horizontal velocity determining the shear rate is the main movement characteristic of the particles, and the vertical velocity indicates the dilatancy phenomenon. Figure 2a shows the time history plot of the equivalent shear strain rate ̇ε e xx . Throughout the evolution process, ̇ε e xx first rapidly increases and then gradually approaches a constant value of 1.0. Figure 2b shows the horizontal velocity distribution contours at six representative moments ( P 1 -P 6 ) selected in Fig. 2a. Overall, under continuous shear, the local flow near the top plate is transferred downward layerby-layer and then gradually expands into the whole system; ultimately, a global steady flow is formed. The horizontal

Horizontal Flow Characteristics
velocity decreases from the top layer to the bottom layer, which is the result of energy dissipation caused by interparticle friction and inelastic collisions. In Fig. 2c, the average horizontal velocity of particles in the same layer is plotted in solid line as a function of the vertical height. Additionally, Fig. 2e shows the time history plots of the average horizontal velocity within some representative layers of particles.
During the evolution process, we divide the whole granular model into flowing and nonflowing regions, written as where H F (t) and H S (t) are the heights of the flowing and non-flowing regions, respectively. Herein, a specific layer is considered flowing only when the average horizontal velocity of particles in this layer is higher than a threshold value, which is set as 1% of the top plate's sliding velocity. We further normalize the height of the flowing region using DH F = H F (t)∕H(t) and present its evolution as a function of ̇ε e xx in Fig. 2d. The evolution results regarding the horizontal flow characteristics in Fig. 2a-e directly reflect the change in flow states. During the time interval from P 1 (0.5 s) to P 3 (6.5 s), the flowing region expands from the local region near the top plate to the global system, and the horizontal velocity of , and the local inertial number (h) along the height the already flowing region gradually increases. Correspondingly, ̇ε e xx increases from 0 to approximately 0.6, and DH F quickly increases from 0 to 1.0. It is worth mentioning that the system at P 3 (6.5 s) has achieved a global flow since DH F = 1 in Fig. 2d. During the time interval from P 3 (6.5 s) to P 6 (27 s), the granular system maintains a global flow state, meanwhile the horizontal velocities gradually increase to constant values for all the selected layers. Moreover, by comparing the evolution characteristics, it is found that the closer to the top plate, the shorter time is required for the particle layer to reach a steady state. In other words, the growth rate of the horizontal velocity decreases from top to bottom, which is caused by the layer-by-layer attenuation of shear momentum. At P 6 (27 s), the system has reached a steady flow state, under which ̇ε e xx approximately reaches 1.0 (Fig. 2a) and the horizontal velocity profile is almost linearly distributed in the vertical direction (Fig. 2c).
Therefore, according to the evolution characteristics of the horizontal flow, the flow status experiences three stages, namely, transmission (T), adjustment (A), and stabilization (S). As outlined in Fig. 2a  In addition, we employ the distribution profile of the local inertial number along the height direction to characterize the flow state. Considering the variation in the shear rate with height, the local inertial number I(y) is calculated as follows: In Fig. 2f, when gravity is absent, the normal pressures along the height P(y) fluctuate around 50 Pa, which is caused by the applied normal stress on the top plate. As a result, the profile shapes of the velocity gradient (Fig. 2g) are similar to those of the local inertial number (Fig. 2h). Interestingly, during evolution, I(y) quickly decreases in the top region but

Vertical Velocity and Shear Dilatancy
Dilatancy is the volume expansion observed in granular materials when they are subjected to shear deformation. Figure 3 illustrates the shear dilatancy mechanism of dense granular flow. Initially, the particles are related to each other through contacts and form a dense packing with the maximum solid fraction, as shown in Fig. 3a. When the solid fraction is larger than a critical value, namely the jamming point, a preliminary dilation is necessary for shear deformation to occur. In other words, for rigid particles, if the system's height undergoing shear remained fixed, the granular system would be blocked. To avoid blockage, the volume has to be expanded. As shown in Fig. 3b, when sheared by the top plate, particles in the flowing region cross the peaks and valleys formed by the neighboring particles, causing the microscopic arrangement of particles. Whether a particle is crossing a peak or valley is indicated by the sign of the particle's vertical velocity. It is noteworthy that a particle could repeatedly cross peaks or fall into valleys during evolution, resulting in vertical velocity fluctuations. Moreover, as the horizontal velocity increases, the duration for a particle horizontally moving over a valley is shortened; hence, the vertical distance falling into the valley is reduced, leading to an increase in dilatancy. Figure 4 shows the vertical velocity distribution contours of all the particles at six different moments in the evolution process. Additionally, the average vertical velocity of particles in the same layer along the height direction is plotted on the right side for each moment to quantify the local dilatancy rate. First, from 0.0001 to 4.6 s, it is observed that the upward and downward moving particles randomly coexist within the system, and the vertical movement gradually propagates from top to bottom. At t = 0.0001s , the average vertical velocity is positive over the entire height, suggesting that the granular system dilates globally. In particular, the particles near the top plate expand rapidly under shear. Soon after that, the global dilatancy is locally replaced as shown at t = 0.1s . At t = 1s and t = 2.7s , the dilating region mainly resides near the top plate ( y∕d > 80 ), as the current dilatancy mainly comes from the local shear flow of particles. At t = 4.6s , the maximum velocity has been transferred from near the top plate ( y∕d > 80 ) to the middle part (y/d = 50-80), which implies that the dilating region is moving down due to the propagation of local flow. Moreover, the sign of the average vertical velocity near the top plate changes from positive to negative. After the granular system enters a global flow state, the profile of vertical velocity along the height approaches zero, and the amount of shear dilatancy reaches the maximum and slightly fluctuates in the subsequent stages.
The dilatancy extent is measured by the normal strain in the vertical direction, yy = H t − H 0 ∕H o , where H 0 is the initial height of the system and H t is the height of the system at time t . Figure 4b and c show the evolution results of yy as a function of time t and of equivalent shear strain rate ̇e xx , respectively. In the evolution process, yy gradually increases to the maximum value of approximately 0.05 and then stabilizes. Moreover, the increase in yy mainly occurs in the transmission stage ( 0s ≤ t ≤ 6.5s, 0 ≤̇e xx ≤ 0.62 ), indicating that shear dilatancy is a prerequisite for particle flow.

Interfacial Friction Coefficient
This part continues to use the simulation results under the condition of P = 50 Pa and V s = 0.5 m∕s to analyze the transient evolution of the boundary driving strength. Figure 5 shows the evolution curves of the scaled normal stress n ∕P and the scaled tangential stress t ∕P , all of which act upon the top plate. To maintain the applied external normal stress constant throughout the evolution process, the position of the top plate must be dynamically adjusted. In detail, once the scaled normal stress n ∕P on the top plate exceeds 1.0, the top plate moves slightly upward to reduce the normal stress, and vice versa. Therefore, as shown in Fig. 5a, n ∕P always fluctuates around 1.0. Additionally, the fluctuation extent is closely related to the flow state, in the way that the fluctuation is intensive in the transmission stage (T) but mild in the adjustment (A), as shown in Fig. 5a. As previously discussed, the particle flow state near the top plate changes dramatically at the beginning of the shear, and the horizontal and vertical velocities increase rapidly. Under this case, the particles collide violently with the top plate, resulting in large fluctuations in the normal and tangential stresses. Subsequently, the particle flow state changes smoothly, so the fluctuations gradually stabilize. In Fig. 5c, the interfacial friction coefficient ′ is calculated based on the results of n and t in Fig. 5a and b. First, it is found that the interfacial friction coefficient ′ gradually decreases to a constant value during evolution. In Fig. 5d, as the shear dilatancy extent increases, the solid fraction decreases; in other words, shear dilatancy makes the system loose. Moreover, mainly decreases in the transmission stage (T) and remains at the minimum value afterward. In the transmission (T) stage, dilatancy loosens the system's spatial geometric constraint, thus weakening the tangential load-bearing capacity and making the interfacial friction coefficient ′ drop rapidly, which can be understood as the strain softening of granular material. However, in the adjustment (A) and stabilization (S) stages, the flow velocities within the particle layers change very slightly, so the shear dilatancy rate is close to zero, and the interfacial friction coefficient ′ also tends to be a constant value.

Stress-Strain Behavior and Effective Friction Coefficient
Furthermore, as shown in Fig. 6, we analyze the evolution principles of the internal stress components and their ratios.
Here, the stress components and their ratios show obvious fluctuations during evolution, which are mainly caused by the microscopic rearrangement of particles. In Fig. 6a, the scaled internal normal stress yy ∕P always fluctuates around 1.0. Moreover, due to the shear flow of the particle system, the internal shear stress xy rapidly increases to the maximum value in the beginning ( ̇e xx ≤ 0.1 ), and then decreases slightly and gradually tends to a constant value, as shown in Fig. 6b. In Fig. 6c, xx ∕ yy is always close to 1.0, since As mentioned earlier, the ratio of the internal shear stress to the internal normal stress is defined as the effective friction coefficient e = xy ∕ yy . Similar to xy , e first rapidly increases to the peak and then gradually drops to a constant value. Comparing Figs. 5 and 6, the changes in n ∕P and yy ∕P are basically the same, both fluctuating around 1.0. However, before reaching a steady state, the internal shear stress xy is always significantly less than t . As a result, in Fig. 6d, there is a significant difference in the evolution trend between the effective friction coefficient e and the interfacial friction coefficient ′ . Moreover, ′ is greater than e until they become approximately equal in the steady flow state. Therefore, the evolution of the flow state from initially static to finally steady under continuous shear, is essentially realized by achieving the final dynamic balance between the boundary driving strength and the system's shear resistance.

Spatial Distribution of the Contact Force Network
The contact force network can indicate the spatial distribution of contact forces and the force transmission paths inside the granular system. Figure 7a clearly presents the contact force network graphs at five different moments under the conditions of P = 50 Pa and V s = 0.5 m∕s . Here, the existence of an interparticle contact is represented by a line segment connecting the center of the two contact particles, and the force magnitude f c is proportional to the line thickness. Initially, the contact force network is dense and uniform. In the evolution process, the contact force network becomes sparsely scattered, and the spatial distribution of strong contacts has a preferred direction. We further conduct a quantitative analysis of the contact force network. During the transient process, N c , N s c , and N w c all gradually decrease with increasing ̇e xx (Fig. 7b-d), which is caused by shear dilatancy. Meanwhile, the coordination number Z slightly decreases from 3.6 to 2.3 during evolution, indicating that the compactness of the filled particles gradually becomes loose. Our data are consistent with the reported Fig. 6 Evolution curves of the scaled internal normal stresses xx ∕P and yy ∕P (a), the scaled internal shear stress xy ∕P (b), the ratio xx ∕ yy (c), and the interfacial friction coefficient ′ and the internal effective friction coefficient e = xy ∕ yy (d), where P = 50 Pa and V s = 0.5 m∕s Page 11 of 18 38 data in reference [23]. Although the numbers of both strong contacts N s c and weak contacts N w c decrease significantly during evolution, the contact number ratios R s c and R w c are basically unchanged at R s c ≈ 0.37 and R w c ≈ 0.63 . In Fig. 7e, the force bearing ratio of strong contacts B s f gradually increases from 0.68 to 0.73 during evolution, and a clear overshoot in B s f occurs at the onset moment of shearing movement. Although the number ratio of weak contacts R w c is always greater than that of the strong contacts R s c , the force bearing ratio of the strong contacts B s f is much greater than that of the weak contacts B w f , indicating that the strong contact force network carries a larger portion of the external load. Therefore, the evolution of the flow state changes the microstructure characteristics of the contact force network and improves the load-bearing capacity.

Anisotropies in Microstructure and Internal Shear Strength
Based on the contact force network in Sect. 5.1, we are ready to deduce the internal shear strength and anisotropy parameters using the stress-force-fabric relationship. First, as shown in Fig. 8a-c, the three distribution functions P( ) , f n ( ) , and f t ( ) at four representative moments are compared to those of the initial state. Initially, P( ) is approximately circular as a c ≅ 0 , indicating that the contact orientations are distributed isotropically. However, the distributions of both normal and tangential contact forces are anisotropic with a n = 0.237 and a t = 0.075 ; particularly, the anisotropy in the normal contact force is very strong. Without flow, the contact forces are mainly used to carry the applied normal stress, while the internal shear strength is relatively weak so that it is easy for the system to flow. In the evolution process, obvious anisotropy is introduced into the distributions of contact orientation and contact forces. Moreover, their preferred directions are collinear, i.e., c ≈ n ≈ t ≈ θ , and here, θ ≈ − 45°. As a result, with the increase in ̇e xx , a c first quickly increases from 0 and then gradually reaches a constant value of approximately 0.32 (Fig. 8a). After passing a peak, a n gradually decreases to a constant value of approximately 0.24 (Fig. 8b). Apart from the initial abrupt drop and rebound, a t basically fluctuates around a small value of 0.06 (Fig. 8c). In the steady state, the ratio of a c , a n , and a t is approximately 16:12:3, implying that the anisotropies in the contact orientation and in the normal contact force are the main sources of the internal shear strength.
Next, the reasons for the changes in the anisotropy parameters are explained in combination with the system's internal micromechanical behaviors. In the beginning of shear, the increase in a c and a n could be attributed to the potential shear resistance caused by the spatial geometric constraints. Nevertheless, most particles except a few near the top plate have not undergone shear flow, so the tangential contact force anisotropy parameter a t changes slightly. Since the horizontal flow is the predominant deformation pattern, the differences in the tangential contact force are slight; thus, a t is low and makes a small contribution to the internal shear strength. Subsequently, because the occurrence of dilatancy reduces the spatial geometric constraints among particles, the original particles that bear large normal contact forces lose lateral support and cannot maintain stability. Once these particles are out of balance, the breakage of the established contacts causes a decrease in a n . On the other hand, particle flow reduces the number of contacts, especially in the direction of the minimum principal stress, which leads to a gradual increase in the contact number difference between along Fig. 8 a-c The transient evolution of anisotropy parameters a c , a n , and a t , as well as the three distribution functions P( ) , f n ( ) , and f t ( ) at representative moments, and (d) * and e as a function of ̇e xx Fig. 9 a Time history plots of the dimensionless equivalent shear strain rate ̇e xx , and b the evolution curves of yy as a function of ̇e xx under four different velocities the direction of the maximum principal stress and along the direction of the minimum principal stress; thus, a c gradually increases. In short, the increase in a c leads to a reduction in a n , which damages the system's stability.
Last, the internal shear strength * is obtained according to Eq. 6, as shown in Fig. 8d. When combining the contributions of the anisotropic parameters, * first rapidly increases to the peak and then gradually drops to a fixed value. Moreover, the transient evolution curves of * and e are basically consistent, which confirms that the internal shear strength originates from the anisotropies in both contact orientation and contact forces.

Effects of Shear Velocity on the Transient Flow State
We select four velocities of the top sliding plate (0.5, 1.0, 1.5, and 2.0 m/s) to study the effects of shear velocity on the rheological properties. The evolution profiles of the horizontal velocity at the other velocities are similar to those of V s = 0.5 m∕s . In addition, three distinct stages, i.e., transmission (T), adjustment (A), and stabilization (S), are identified for all velocities. Nevertheless, there still exist differences in the evolution of the flow state among them. As shown in Fig. 9a, with increasing shear velocity, the time required for ̇e xx to reach the maximum value increases, indicating that the evolution rate of the flow state decreases. Under inertial flow, an increase in the shear velocity strengthens the inertial effect and reduces the duration of contact between particles, thereby weakening the momentum transfer between particles. Therefore, the time required for reaching steady state becomes longer with increasing shear velocity. Additionally, in Fig. 9b, the maximum normal strain yy , i.e., the dilatancy extent, increases with shear velocity. Therefore, the shear velocity does not affect the evolution characteristics of particle flow, but changes the momentum transfer rate.

Effects of Shear Velocity on Macroscopic and Microscopic Mechanical Properties
We present the evolution results of normal stress and shear stress acting on the top plate ( n and t in Fig. 10a and b) and within the granular system ( yy and xy in Fig. 10c and  d), respectively. In terms of the evolution characteristics, these stresses fluctuate more intensively as the shear velocity increases, especially for the boundary stresses n and t . In addition, Fig. 10e shows that the internal shear stress  xy significantly increases with increasing ̇e xx at higher shear velocities. On the other hand, with increasing shear velocity, the magnitudes of shear stresses t and xy at the same moment significantly increase, while the scaled normal stresses n ∕P and yy ∕P always fluctuate around 1.0. As a result, both the interfacial friction coefficients ′ and the effective friction coefficients e increase as the shear velocity increases, as shown in Fig. 10c and f. This could be explained by the fact that when the overall shear rate increases, the boundary driving strength must be correspondingly improved to overcome more energy dissipation and maintain steady flow.
Furthermore, the anisotropy parameters a c ,a n , and a t all increase with increasing shear velocity as shown in Fig. 11a-c. According to Eq. (6), the system's internal shear strength characterized by * is correspondingly improved in Fig. 11d. The changes in these anisotropy parameters could be explained by comparing the contact force network graphs in the steady state. As shown in Fig. 12a, at V s = 0.5 m∕s , a dense and continuous contact force network is formed as the particles remain in contact and interact frictionally with neighbors over long durations. By contrast, at V s = 2 m∕s , the distribution of contact force is so dispersed that only some local line segments of contact force exist. As a result, the difference in contact number increases between along the maximum principal stress and along the minimum principal stress, leading to an increase in the contact orientation anisotropy parameter a c . As shown in Fig. 12b, in the steady state, both the strong contact number N s c and the weak contact number N w c decrease with increasing shear velocity. N w c is always greater than N s c . Additionally, with increasing shear velocity, the bearing ratio of weak contacts B w f decreases while the bearing ratio of strong contacts B s f increases, so the differences in the force magnitudes become more distinct, which explains the increase in the anisotropy parameters of contact forces a n and a t .
Finally, we present and compare the different types of friction coefficients in the steady state, as shown in Fig. 13a. Numerous experimental and simulation results have confirmed that for the inertial regime, the macroscopic friction coefficient increases with an increasing inertial number I in the following manner: where s , g , and I 0 are material-dependent constants. In Fig. 13a, our simulation results can be well fitted using Eq. (8), where the fitting coefficients and confidence are listed in Table 2. These fitting equations form the constitutive laws for the dense granular material under plane shear in this study. The three types of friction coefficients are approximately equal because the boundary driving strength and the system's shear resistance are dynamically balanced.
Furthermore, Fig. 13b shows the variations in the anisotropy parameters with the inertial number I. Specifically, the contact orientation anisotropy parameter a c increases significantly from 0.34 to 0.47, and the normal contact force anisotropy parameter a n gradually increases from 0.23 to 0.30, while the tangential contact force anisotropy parameter a t only slightly increases from 0.06 to 0.09. By comparing their increasing extents, it can be concluded that a c and a n are the main sources of the system's shear strength. Moreover, the contribution of a c becomes more predominant with the increase in the inertial number I.  Fig. 13 a The interfacial friction coefficient ′ , the effective friction coefficient e and the internal shear strength * , and b the variations in the anisotropy parameter a c , a n , and a t , as a function of the inertial number I Table 2 Fitting coefficients of ′ , e and * using the form of = s + ( g − s )∕(I 0 ∕I + 1)

Conclusions
DEM simulations are employed to investigate the transient rheological properties of a dense granular system under plane shear in the inertial regime. Simulation results show that there are three typical stages in the evolution process of the system's horizontal flow state, namely, transmission, adjustment, and stabilization. The shear dilatancy which is caused by the vertical movement of particles, weakens the system's tangential load-bearing capacity, decreases the interfacial friction coefficient ′ and reduces the boundary driving strength. On the other hand, the shear flow induces variations in the anisotropies of the contact orientation and contact forces, thus improving the system's internal shear strength * and increasing the effective friction coefficient e . As a result, the evolution trends are opposite between ′ and e ; moreover, ′ is greater than e until approximately equal in the steady flow state. Therefore, the flow evolution from static to steady flow is realized by achieving the final dynamic balance between the boundary driving strength and the system's shear resistance. In addition, an increase in the shear velocity enhances the rheological evolution rate, but the evolution pattern is not changed. Distinguished from the micromechanical behaviors, the system's shear strength always mainly originates from the anisotropies in contact orientation and in normal contact force. Moreover, the anisotropy in contact orientation makes more contribution as the shear velocity increases.
Our work comprehensively analyzes the competition between the boundary driving strength and the system's shear resistance in the evolution process, and thus reveals the shear equilibrium mechanism of dense particle flow. In addition, the effects of shear velocity on the transient evolution are determined, by means of integrated analyses on both the macroscopic quantities and the microstructure characteristics. Together with the previous studies on steady flow, our transient analysis would help supplement a more complete understanding of the flow behaviors of granular materials. In the future, we plan to distinguish the contributions of strong contacts and weak contacts on the granular flow behaviors.