HTML
-
To simulate the dynamical behaviors of millions of particles, this work was carried out on multiple GPUs by the Discrete Element Method (DEM) code we developed[16-17]. In the DEM model the interactions between the particles are given by Hertz-Mindlin contact model[18-19]. Assuming the particles are identical, the normal and tangential contact forces between two contacting particles are:
where G is the shear modulus, E is the Young’s Modulus, r is the radius of the particles and m is the mass of particles.
$ {{\delta }}_{{ij}_{\rm{n}}} $ and$ {{\delta }}_{{ij}_{\rm{t}}} $ are normal and tangential displacement vectors, and$ {\delta }_{{ij}_{\rm{n}}} $ and$ {\delta }_{{ij}_{\rm{t}}} $ are their modules, respectively.$ {{v}}_{{ij}_{\rm{n}}} $ and$ {{v}}_{{ij}_{\rm{t}}} $ are normal and tangential relative velocities between the particles I and j.$ \kappa $ is related to the coefficient of restitution. If there is friction, the Coulomb yield criterion |${ F}_{{ij}_{\rm t}} $ | ≤ |μs${ F}_{{ij}_{\rm n}} $ | is satisfied by truncating the magnitude of${ F}_{{ij}_{\rm t}} $ . As a result, if |${ F}_{{ij}_{\rm t}} $ | > μs|${ F}_{{ij}_{\rm n}} $ |,${ F}_{{ij}_{\rm t}} $ =μs|${ F}_{{ij}_{\rm n}} $ |uij/|uij|[19].${ F}_{{ij}_{\rm t}} $ is the friction between two contacting particles i and j, and μs is the friction coefficient.Under gravity field, the equations of motion of the particles are
These equations are solved by integration using the Velocity-Verlet scheme[20].
-
To simplify the target model, here we investigated the granular flow down an inclined chute discharging through the opening gate. A system, consisting of 1 000 000 mono-disperse beryllium spherical particles (diameter d is 5 mm) was used and the material properties of the reservoir and chute was the same as particles. The walls of the reservoir and chute were smooth and frictional. A Cartesian coordinate system was established: the flow direction was set to x axis and x=0 was at the opening gate; the span-wise direction was set to y axis and y=0 was along the center of the chute; the height (z) direction was perpendicular to the chute and z=0 was at the bottom of the chute. The configuration of the model is schematically shown in Fig. 1 and the material parameters are presented in Table 1. The time-step in DEM simulations was 5×10–7 s. Initially a random filling process was used to fill the reservoir with the particles until reach a static packing. Then the gate was opened and steady state flow was allowed to develop before acquiring the data used for subsequent analysis.
Quantity Symbol/Unit Value Diameter of particles d/mm 5 Width of the chute W/mm 40d Length of the chute L/mm 600d Height of the gate G/mm 25d Inclination angle of the chute θc/(°) 25° Length/width of the reservoir D/mm 40d Elastic modulus E/GPa 287 Poisson’s ratio ν 0.032 Particle-particle and particle-wall friction coefficients μ 0.2 Density ρ/(kg/m3) 1 850 Particle-particle and particle-wall coefficient of restitution ε 0.9 Hardness H/GPa 1.67 Table 1. Geometrical and material parameters.
-
Here we use granular flows down inclined flat planes, with periodic boundaries in the x (flow) and y (span-wise) directions, to study the stability of the inclined flows. A system, consisting of 10 000 mono-disperse beryllium spherical particles (diameter d is 5 mm) was used. A Cartesian coordinate system was established: the flow direction was set to x axis and x=0 was at the center of the length; the span-wise direction was set to y axis and y=0 was along the center of the chute; the height (z) direction was perpendicular to the chute and z=0 was at the bottom of the chute. The parameters of material properties are also shown in Table 1. The time-step in DEM simulations was 5×10–7 s too. The simulation cell dimensions in the x and y direction were 20d and 20d respectively, resulting in an approximately 25d height in the z direction. Initially, 10 000 particles randomly dropped under gravity into the horizontal plane and a static packing was achieved after a while. Then the system was tilted to a desired inclination angle (25°) to initiate the flow. At this inclined Angle the granular flow will be accelerative. After reaching a steady flows, perturbation was added to test the stability of the flows and then Lyapunov exponent
$ \xi $ was calculated by fitting the formula:$\Delta E\left(t \right) = \Delta E\left(0 \right){\rm{exp}}\left({ - \xi t} \right)$ , where$ \Delta E $ is the absolute average difference of per particle between the perturbed system and the control no-perturbed system. In our study, three perturbation protocols (Fig. 2) were considered: (a) changing the gravity from 1 g to 5 g during t=0~0.1 s; (b) changing the gravity from 1 to –1 g during t=0~0.1 s; (c) pining the bottom particles instantaneously at t=0 s. -
In the process, Voronoi volume of each particle was calculated and the local volume fraction was ratio of particle volume and its Voronoi volume. Then these local volume fractions were aligned as property of the corresponding particle and it was time-averaged by the same procedure of the particle velocity[21]. In the procedure, hundreds of snapshots were gathered and the space of each snapshot was divided into the same meshes. The particle quantities (including velocity, volume fractions and etc.) of each mesh were averaged by the appearance time of particles within this mesh to obtain the time-averaged quantities.
The virial stress is commonly used to find the macroscopic (continuum) stress in molecular dynamics systems, and is defined as
where
$ k $ and$ l $ are atoms in the domain,$ V $ is the volume of the domain,$ {m}^{\left(k\right)} $ is the mass of atom$ k $ ,$ {{u}_{i}}^{\left(k\right)} $ is the$ {i}^{th} $ component of the velocity of atom$ k $ ,$ {\bar{u}}_{j} $ is the$ {j}^{th} $ component of the average velocity of atoms in the volume,$ {x}_{i}^{\left(k\right)} $ is the$ {i}^{th} $ component of the position of atom$ k $ , and$ {f}_{i}^{\left(kl\right)} $ is the$ {i}^{th} $ component of the force applied on atom$ k $ by atom$ l $ . In this paper, the time-series profiles of particle positions and velocities were collected and space in the chute was divided into subdomain to calculate the time-averaged virial stresses.
2.1. Contact model
2.2. Simplified model of inclined chute flows
2.3. Perturb the inclined flows
2.4. Data processing
-
In an inclined chute with a rough base, the particles behave like solid, fluid or gas by varying the inclination angles[22-23]. If the base is sufficiently smooth, the flow is simple and adjustable. To simplify the model, monosize-distribution tungsten particles (diameter d is 5 mm) were simulated. The material of the reservoir and chute was the same as particles and their roughness was assumed as nil. The configuration of the ‘the basic’ case is schematically shown in Fig. 3. The geometrical and material parameters in the basic case are shown in Table 1. The time-step in DEM simulations was 5×10–7 s too. The inclination angle was set to 25° to avoid potential stagnant or gaseous state in the downstream chute[23]. Initially the particles were packed in the reservoir and then were discharged through an opening gate. After 20 s later, the time-averaged velocity and volume fraction fields along the chute flow were plotted and the stability of the flow was analyzed.
Figure 3. (color online)Configuration of the simple model of chute flows. The red arrow denotes the flowing direction.
The morphology of the flow was nearly stable over time [see in Fig. 4(a)]. From the gate to x=600d, fluctuation of the height of the accelerating flow was less than 10% which was consistent with Holyoake et al[24]’s results (see in appendix for detailed analysis of the fluctuation and stability).
Figure 4. (color online)The simulation results of granular flows in a chute with gate height of 25d, chute width of 40d and inclination angle of 25°.
The time-averaged velocity and volume fraction fields are shown in Fig. 4(b) and 4(c), where the volume fraction fields were calculated by plotting the Voronoi diagrams[25]. The volume fraction in the chute flow varied from 0.5 to 0.6 except at the free surface and bottom, which was consistent with the volume fraction in random packing. Fig. 4(d) shows the influence of the walls on the volume fraction. Similar to the bottom, the volume fraction near the walls varied between 0.4 and 0.5. From the spatial profiles of
$ {v}_{x} $ , it appeared that there is nearly a plug flow in the chute as the variation of$ {v}_{x} $ is small along the depth. The influence of the walls on$ {v}_{x} $ is shown in Fig. 5(d) and it was found that central particles move more quickly than particles near the sidewalls (the discrepancy is less than 10%).Figure 5. (color online)Spatial profile of the simulation results with gate height of 25d, chute width of 40d and inclination angle of 25°.
From former studies, it showed there is a basal rolling layer in a stable, fully developed flow on the flat, frictional incline[26], which was also observed in this accelerating flow [see in Fig. 5(c)]. Interestingly, the result showed that the wy of second layer is opposite, which agrees with the previous picture[27]. The profile of width-averaged and depth-averaged vx could be fitted by using the following function:
in Ref. [24], fitting parameters
$ {v}_{0}=1.191\;{\rm{m}}/{\rm{s}} $ ,$ \beta =4.488\;{\rm{m}}/{{\rm{s}}}^{2} $ ,$ \gamma =0.024\;{{\rm{m}}}^{-1} $ , which implied that vx will converge to 13.7 m/s when$ x\to \infty $ . The height h(x) could also be fitted if the flux conversation was assumed in the accelerating flow:Here the fitting parameter
$ {C}_{0}=0.113\;{{\rm{m}}}^{2}/{\rm{s}} $ and fitting result is shown in Fig. 4(a). The influence of sidewalls on the velocity profile could be evaluated by using fitting formula in Ref. [28]:It is noted that this formula is used for studying about the influence of sidewalls on the surface velocity in granular avalanche[28] but it is also available for inner and bottom velocities here [see in Fig. 5(d)]. The fitting parameter
$ \varLambda =0.037\;{\rm{m}} $ here and was an invariant. To evaluate the stability of the flow, the temporal profiles of velocity and volume fraction fields were analyzed and their standard error was calculated. Here it was showed that the standard errors of volume fraction are very small except at the free surface. For$ {v}_{x} $ , the standard error on the bottom was more than the standard error at the free surface and both were relatively small.To study the inferences of material parameters on the chute flows, comparing simulations were performed and the results are shown in Fig. 6 and Table 2. It appeared that the friction coefficient and chute width influence the flows. The flow rate was bigger when the particles had a low modulus of elasticity (E=5×106 Pa), which was consistent with previous observation in hopper flow[29]. Moreover, the flow was slower in the narrower chute.
Case Mass flow rate/(kg/s) Basic case 26.8 Softer particles (E=5×106 Pa) 29.3 Higher friction coefficient (μ=0.5) 18.8 Narrower chute (W=10 cm) 10.6 Bigger particles (d=10 mm) 25.6 Higher density of particles (ρ=7 800 kg/m3) 120.4 Table 2. The influences of material and geometrical parameters on the mass flow rate.