# Particle Directional Conveyance under Longitudinal Vibration by considering the Trough Surface Texture: Numerical Simulation Based on the Discrete Element Method

Authors: Huazhi Chen; Shengyuan Jiang; Rongkai Liu; Weiwei Zhang
Journal: Shock and Vibration
 (2018)

---

## Abstract

Particles can move directionally in a trough with finlike asperities under longitudinal vibrations. Here, we present an analysis of the particle conveyance mechanism and the influence of the asperity shape on the particle conveyance capacity by employing a numerical simulation based on the discrete element method (DEM). A dynamic-static matching method is proposed to characterize the three microcontact parameters in the simulation: the restitution coefficient, static friction coefficient, and rolling friction coefficient. The simulation shows that the asymmetric force induced by the finlike asperities and its cumulative effect over time lead to the particle directional conveyance. The conveyance velocity increases with increasing vibration time and is related to the median coordination number. The asperity height and slope inclination angles determine the trough shape and distance between two asperities directly. An undersized or oversized distance reduces the steady conveyance velocity. We find the optimal distance to be between one and two particle diameters.


---

## Body

## 1. Introduction

In situ resource utilization (ISRU) is an important basis for establishing a future lunar base. ISRU involves using lunar regolith and rock as building materials and manufacturing metal and oxygen from these raw materials [1–3]. In this process, conveyance is essential for raw material, particularly granular materials (e.g., lunar regolith), as shown in Figure 1. In most cases, the granular materials must be transported from the ground to a certain height automatically and continuously. Considering the constraints of lunar engineering, a vibratory conveyor with a simple configuration that consumes relatively little energy is more suitable than a belt conveyor [4]. Furthermore, a vibratory conveyor can transport high-temperature or corrosive materials because its trough can be readily made from a refractory and corrosion-resistant material. These properties make such conveyors widely preferred in the metallurgical and chemical industries. In view of the above advantages, several investigations have been conducted to extend vibratory conveyance theory and improve the performance of conveyors in engineering. Two types of vibratory conveyor are currently applied in practice and are classified according to the vibration direction relative to the trough surface. The first features sinusoidal vibration acting obliquely on the trough [5–8], and the second features vibration parallel to the axis of the trough (i.e., the longitudinal direction) but with a nonsinusoidal profile [9, 10]. Materials cannot be transported directionally if the trough vibrates longitudinally in a sinusoidal form. However, such transport is enabled if the trough has an appropriate surface morphology, for example, a sawtooth texture [11]. We previously reported such behavior [12]. This topic merits comprehensive study because of its potential value for improving the performance of current conveyors.Figure 1
Conception of lunar regolith excavation and conveyance to establish a lunar base [2].Generally, granular materials are considered to be a point mass for convenient theoretical analysis, but this simplification cannot reflect the actual properties of such materials. Previous research has established that granular materials exhibit complex properties under vibration, such as convection, heaping, size separation, and even climbing in narrow channels [13–16]. A large number of particle interactions leads to a kinestate distinct from that of a single body on a vibratory conveyor. The presence of a variety of multiple types of particles, particularly tiny and opaque particles, makes it challenging to experimentally study detailed particle motion in a trough. Fortunately, numerical simulations based on the discrete element method (DEM), introduced by Cundall and Strack [17], can be used to study this configuration. The advantage of the DEM is that it provides detailed particle-scale information for analysis. DEM simulations have provided an effective means to study granular material behavior over several decades of development [18–22].In our previous work, we studied particle directional conveyance in a trough with sawtooth asperities using DEM simulations [12]. The primarily goal was to study the influence of the vibration conditions on the particle conveyance capacity. For the DEM simulations, the accuracy of the parameters set for the calculation can significantly affect the reliability of the result. In the simulation, we characterized particle microcontact parameters such as the static friction coefficient by matching the angle of repose (AOR) [23, 24] with the experimental value. This method is reliable when considering particles with static or quasi-static characteristics but is not completely suitable for studying the dynamic characteristics, which is an unavoidable issue as long as we employ DEM simulation in this context.For the trough surface morphology, Farkas et al. [25] observed that the size of the finlike asperities in an annular channel can affect particle horizontal conveyance under vertical vibrations. Thus, we chose to explore the influence under a longitudinal vibration. Distinct from our previous work, in this study, we focused on the mechanism of particle directional conveyance and the influence of the finlike asperity size on the conveyance capacity using DEM simulations. To improve the reliability of the simulation results, we propose a dynamic-static matching method to determine the particle microcontact parameters. This paper is organized as follows. In Section 2, we introduce the DEM model used to calculate the contact force and torque. In Section 3, we review the simulation conditions and the dynamic-static matching method to characterize the particle microcontact parameters. In Section 4, we present and discuss the simulation results. Finally, we report our conclusions in Section 5.


## 2. DEM Model Description

In our simulation, a three-dimensional DEM based on the soft-sphere model [26–28] was used. The particles were considered to be dry and clean; that is, there were no cohesive forces and no liquid bridges between the particles. Each particle exhibited two types of movement, translational motion and rotational motion, as governed by Newton’s second law:(1)midvidt=mig+∑j=1niFij,Iidωidt=∑j=1niTij,where vi and ωi are the translational and angular velocity vectors of particle i, respectively; mi and Ii are the mass and the rotational inertia of particle i, respectively; g is the gravitational acceleration; ni is the number of particles j in collision with particle i; and Fij and Tij are the contact force and torque generated between particles i and j, respectively, and are composed of several components: namely,(2)Fij=Fn,ij+Ft,ij,Tij=Tn,ij+Tt,ij+Tr,where Fn,ij and Ft,ij are the normal and the tangential forces, respectively, and Tn,ij and Tt,ij are the torques caused by the normal force and the tangential force, respectively. Tr is the rolling friction torque when relative rotation occurs between the contacting particles.Due to the absence of cohesive forces and liquid bridges, we chose the no-slip Hertz-Mindlin contact model [29] to calculate Fij and Tij. This model is based on the classical Hertz theory along the normal direction [30] and the work of Mindlin [31] along the tangential direction. Depending on the model, the contact forces that include particle-particle and particle-wall interactions are as follows:Normal force(3)Fn,ij=-43E∗R∗δn,ij3-256βcSnm∗vn,ijnij.Tangential force(4a)Ft,ij=-Fte,ijn-256βcStm∗vt,ijforFt,ij≤μsFn,ij-μsFn,ijvt,ijvt,ijforFt,ij>μsFn,ij,(4b)Fte,ijn=Fte,ijn-1+StnΔδt,ijforΔFn,ij≥0,(4c)Fte,ijn=Fte,ijn-1StnStn-1+StnΔδt,ijforΔFn,ij<0.Torque(5)Tn,ij=Rinij×Fn,ij,Tt,ij=Rinij×Ft,ij,Tr=-μrFn,ijRiωni.The parametersμs and μr are the static and rolling friction coefficients, respectively; δn,ij and δt,ij are the normal overlap and tangential displacements, respectively; vn,ij and vt,ij are the relative normal and tangential velocities, respectively; nij is the unit vector pointing from the center of particle i to particle j; ωni is the unit angular velocity vector when a relative rotation occurs; and E∗, R∗, and m∗ are the equivalent Young’s modulus, particle radius, and mass, respectively [32].In (4a), (4b), and (4c), the superscript n is the time step number, Fte,ij is the elastic component of the tangential force, and Δδt,ij and ΔFn,ij are the tangential displacement increment and the normal force increment, respectively. Fte,ij depends on the interaction phase. If the interaction is in the loading phase, the increment of Fte,ij is simply proportional to Δδt,ij (see (4b)). If the interaction is in the unloading phase, then the elastic tangential force at the previous time step should be reduced (see (4c)) due to the reduction in the contact area during unloading. This modification is intended to prevent the creation of spurious energy [33].The parametersSn, St, and βc are given by(6)Sn=2E∗R∗δn,ij,St=8G∗R∗δn,ij,βc=ln⁡eln2e+π2,where G∗ is the equivalent shear modulus and e is the restitution coefficient.Based on the no-slip Hertz-Mindlin contact model and the matched calculation procedure, the particle’s kinestate in each time step is obtained. In this work, the EDEM software package [34] was employed for the numerical simulation. However, before calculation, the conditions should be set in detail, particularly the microcontact parameters e, μs, and μr, as described in Section 3.


## 3. Simulation Setting

### 3.1. Three-Dimensional Particle Model

The particle was defined as a glass bead of 0.4 mm to 0.6 mm (average 0.5 mm) diameter and 2,468 kg/m3 real density. This bead presented a nearly perfect spherical shape according to the micrograph in Figure 2 (a reference rod was used to estimate the glass bead diameter). Thus, we employed a spherical particle with a uniform diameter of 0.5 mm for the simulation, and we tentatively ignored the effect of particle size variation in this study.Figure 2
Micrograph of a glass bead to determine the particle model for the simulation.


### 3.2. Microcontact Parameter Matching

The primary challenge in DEM simulation is to determine the material mechanical parameters and the microcontact parameters. The material mechanical parameters, which include the density, Young’s modulus, and Poisson’s ratio, can be determined experimentally or by consulting references. However, the microcontact parameters cannot be obtained accurately using this method. Generally, such parameters are determined by matching a parameter that is measured easily in both simulation and experiment [35, 36]. In this study, the objective parameter we used was the conveyance rate (Rc), which is defined as the ratio of the transported volume per second to the total volume. Relative to matching the stress in a triaxial compression test [37] or the AOR in a sandglass test [24], this method is better suited to study the particle dynamic motion and obtain more accurate microcontact parameters. The main procedure of the method is described in Figure 3, and the details are introduced in Figure 3.Figure 3
Procedure to determine the microcontact parameters for the DEM simulation.Figure4 shows a sketch of the experiment conditions. The trough was made from synthetic resin and augmented with finlike asperities along a surface measuring 100 mm × 10 mm (LT×WT), as shown in Figure 4(a). The trough was positioned obliquely to obtain an inclination angle (θ) of 10°. The height of each asperity (h) was 0.3 mm. The tangent values of the slope inclination angles α and β were 6/17 and 2, respectively. A silo with an inner diameter (d0) of 3 mm was placed under the discharge hole to evaluate the conveyance capacity. In addition, a cover plate was attached to the trough to prevent the glass beads from being ejected from the trough during vibration.Figure 4
Sketch of the experimental conditions. (a) The trough with finlike asperities and the main experimental process. (b) Experimental equipment for vibratory conveyance.
(a) (b)Vibration was obtained using a vibrator (HEV-50, manufactured by Foneng Technology and Industrial (Nanjing) Co., Ltd., Nanjing, China) with an adjustable amplitude (A) and frequency (f) using a control device, as depicted in Figure 4(b). In this study, the vibration was sinusoidal:(7)x=Asin⁡2πft,where x is the displacement of the trough. The real-time vibration was monitored using a laser displacement sensor (Panasonic HG-C1030).In the experiment, a known volume of glass beads (Vtot = 170 mm3) was placed on the middle of the trough surface containing the finlike asperities (Step 1). Then, the trough was vibrated along its central axis, and the glass beads were collected by the silo when the beads were discharged from the discharge hole (+x-direction in Figure 4, Step 2). The collection process was recorded using a digital camera to ascertain the heaping height (hm). The mass of the glass beads was sufficiently small that its influence on the vibration amplitude A was ignored.WhenA=0.3 mm and f = 50 Hz, the heaping height increased over time until all of the glass beads were discharged, as shown in Figure 5. To evaluate the conveyance rate, we selected an interval during which the relationship between hm and t was approximately linear. We considered that the relationship would be linear if the correlation coefficient (r) was greater than 0.9 in this study. The corresponding conveyance rate is(8)Rc,e=πd024Vtot·dhmdt,t∈6,12.Figure 5
Relationship between the experimental heaping height and time.According to the fitting curve in Figure5, Rc,e was 0.089 s−1 in the experiment.In the DEM simulation, an assembly of 1,500 spherical particles (the total occupation was approximatelyVtot = 170 mm3) was employed and completely generated from the particle generation region in 1 s, as shown in Figure 6. The region measured 20 mm × 10 mm. After generation, the trough vibrated along the x-direction for 10 s. The size of the trough and finlike asperity, the inclination angle (θ), and the vibration parameters (A and f) were the same as in the experiment. To shorten the actual time for simulation, we used the particle number instead of the heaping height to define Rc,s:(9)Rc,s=1Ntot·dNcdt,where Ntot = 1,500 is the total particle number and Nc is the discharged particle number. Similarly, Rc,s was also obtained according to the linear part of the t-Nc curve.Figure 6
Sketch of the DEM simulation conditions.The trough and the particle material mechanical parameters are listed in Table1. The particle Poisson’s ratio was based on Tatemoto et al.’s work [38]. The shear modulus was set to 5×107 Pa to improve the computational efficiency. Though this value is small for glass, the influence can be revised using microcontact parameter matching as described below. Other parameters were obtained from the raw material suppliers.Table 1
Material mechanical parameters of the trough and particles.
Parameter Trough Particles Densityρ [kg/m3] 1070 2468 Shear modulusG [Pa] 8.63 × 10 8 5 × 10 7 Poisson’s ratioγ [-] 0.39 0.23In the DEM simulation, the time stepΔt depends primarily on the particle material mechanical parameters [39]. In this study, Δt was 1.51×10-6 s, and the data were recorded every 0.001 s.It is necessary to analyze the significance of the microcontact parameters before matching. To this end, an orthogonal test is simple and appropriate. Generally, the restitution coefficiente and static friction coefficient μs each lies between 0 and 1, and the rolling friction coefficient μr lies between 0 and 0.1. According to Tatemoto et al.’s work [38], the particle-particle restitution coefficient ep-p is 0.9. Assuming that the static and the rolling friction coefficients of the particle-particle and particle-trough interactions are the same (i.e., μs,p-p=μs,p-t=μs and μr,p-p=μr,p-t=μr), the number of the undetermined microcontact parameters is reduced from 5 to 3. The resulting levels of each parameter in the orthogonal test are listed in Table 2.Table 2
Microcontact parameter levels in the orthogonal test.
Level number (A)ep-t [-] (B)μs [-] (C)μr [-] 1 0.3 0.5 0.08 2 0.5 0.3 0.06 3 0.7 0.9 0.04 4 0.9 0.7 0.10The orthogonal test included 3 factors and 4 levels; thus, 16 subsimulations were required. The relationship betweenNc and t in each subsimulation is shown in Figure 7, and the conveyance rates Rc,s are listed in Table 3.Table 3
Orthogonal test scheme and results.
Subsimulation number (A) Level number (B) Level number (C) Level number R c , s (s−1) 1 1 1 1 0.124 2 1 2 2 0.231 3 1 3 3 0.074 4 1 4 4 0.077 5 2 2 4 0.306 6 2 1 3 0.194 7 2 4 2 0.153 8 2 3 1 0.133 9 3 3 2 0.193 10 3 4 1 0.190 11 3 1 4 0.221 12 3 2 3 0.252 13 4 4 3 0.028 14 4 3 4 0.053 15 4 2 1 0.049 16 4 1 2 0.041  RangeRs 0.684 0.390 0.161Figure 7
Relationship between the discharged particle number and vibration time in each subsimulation: (a)ep-t = 0.3, (b) ep-t = 0.5, (c) ep-t = 0.7, and (d) ep-t = 0.9.
(a) (b) (c) (d)The intuitive analysis (range analysis) in the orthogonal tests [40] is a simple method to estimate the significance of the factors. A greater range (Rs) implies that the corresponding factor is more significant. Based on the result in Table 3, we determined that the restitution coefficient of the particle-trough (ep-t) was the most significant factor and that the rolling friction coefficient μr produced less of an effect on the conveyance rate. In addition, the results of simulations numbers 3 and 4 were more similar to Rc,e. Thus, the parameter ep-t was determined to be 0.3 in this study.Although the influence of the rolling friction coefficient on the conveyance rate was not prominent, this influence should be determined as accurately as possible. An effective method to perform this evaluation is the static matching method, also known as AOR matching in this study. The procedure of measuring the AOR is as follows: (a) place a hopper on the horizontal plane; (b) fill the hopper with a certain volume of particles; (c) hoist the hopper slowly until the particles are completely discharged; and (d) measure the AORθr. Figure 8 illustrates the procedure.Figure 8
Procedure to measure the AOR.The experimental value of AOR (θr,e) was 25.4°, while the simulation values (θr,s) based on simulations numbers 3 and 4 were 28.3° and 34.8°, respectively, as shown in Figure 9. The simulated value approached θr,e when μr=0.04 despite μs being larger.Figure 9
AOR: (a)θr,e = 25.4° in the experiment; (b) θr,s = 28.3° in the simulation (number 3); and (c) θr,s = 34.8° in the simulation (number 4).
(a) (b) (c)To determine the final values ofμs and μr, we employed an optimal test whose object was to match (Rc,s,θr,s) with (Rc,e,θr,e), with the conditions(10)Rc,s-Rc,eRc,e≤ξRc∩θr,s-θr,eθr,e≤ξθr,where ξ(Rc) and ξ(θr) were the allowable errors, both set as 0.05. The results are shown in Figure 10, where μr = [0.03, 0.04] and μs = [0.3, 0.5, 0.7, 0.9]. The figure indicates that Rc,s declined as μs rose, while θr,s increased under the same μr. Based on the curves and allowable errors, the optimized intervals of μs (Aop) under various μr were determined. Clearly, Rc was larger when μr = 0.03, and the ultimate values of μs and μr were determined to be 0.8 and 0.03, respectively.Figure 10
Relationship betweenRc,s, θr,s, and μs under various values of μr.By settingep-p=0.9, ep-t=0.3, μs=0.8, and μr=0.03, the conveyance rate Rc,s and the AOR θr,s were 0.091 s−1 and 26.1°, respectively, which further validated the reliability of these predictions. Then, the total microcontact parameters for the DEM simulation were determined. The entire procedure is termed the dynamic-static matching method.


### 3.3. Dimensions of the Finlike Asperity

The parameters to determine the finlike asperity areh, α, and β, as shown in Figure 4. To simplify our exposition, we introduce the following three dimensionless parameters:(11)Yh=hdp,Ya=tan⁡α,Yb=tan⁡β,where dp is the particle diameter (dp=0.5 mm). The levels of each parameter are listed in Table 4. The scope of the present study did not allow simulation of all possible combinations.Table 4
Dimensionless parameters to describe the finlike asperity and associated levels.
Level number Y h Y a Y b 1 0.2 0.1 1 2 0.4 0.2 2 3 0.6 0.3 3 4 0.8 0.4 4 5 1.0 0.5 5Another reason to use the DEM simulation instead of an experiment in this study is that smaller asperities (particularly whenYh < 0.6 and Ya > 0.3) are challenging to manufacture using conventional methods.


## 3.1. Three-Dimensional Particle Model

The particle was defined as a glass bead of 0.4 mm to 0.6 mm (average 0.5 mm) diameter and 2,468 kg/m3 real density. This bead presented a nearly perfect spherical shape according to the micrograph in Figure 2 (a reference rod was used to estimate the glass bead diameter). Thus, we employed a spherical particle with a uniform diameter of 0.5 mm for the simulation, and we tentatively ignored the effect of particle size variation in this study.Figure 2
Micrograph of a glass bead to determine the particle model for the simulation.


## 3.2. Microcontact Parameter Matching

The primary challenge in DEM simulation is to determine the material mechanical parameters and the microcontact parameters. The material mechanical parameters, which include the density, Young’s modulus, and Poisson’s ratio, can be determined experimentally or by consulting references. However, the microcontact parameters cannot be obtained accurately using this method. Generally, such parameters are determined by matching a parameter that is measured easily in both simulation and experiment [35, 36]. In this study, the objective parameter we used was the conveyance rate (Rc), which is defined as the ratio of the transported volume per second to the total volume. Relative to matching the stress in a triaxial compression test [37] or the AOR in a sandglass test [24], this method is better suited to study the particle dynamic motion and obtain more accurate microcontact parameters. The main procedure of the method is described in Figure 3, and the details are introduced in Figure 3.Figure 3
Procedure to determine the microcontact parameters for the DEM simulation.Figure4 shows a sketch of the experiment conditions. The trough was made from synthetic resin and augmented with finlike asperities along a surface measuring 100 mm × 10 mm (LT×WT), as shown in Figure 4(a). The trough was positioned obliquely to obtain an inclination angle (θ) of 10°. The height of each asperity (h) was 0.3 mm. The tangent values of the slope inclination angles α and β were 6/17 and 2, respectively. A silo with an inner diameter (d0) of 3 mm was placed under the discharge hole to evaluate the conveyance capacity. In addition, a cover plate was attached to the trough to prevent the glass beads from being ejected from the trough during vibration.Figure 4
Sketch of the experimental conditions. (a) The trough with finlike asperities and the main experimental process. (b) Experimental equipment for vibratory conveyance.
(a) (b)Vibration was obtained using a vibrator (HEV-50, manufactured by Foneng Technology and Industrial (Nanjing) Co., Ltd., Nanjing, China) with an adjustable amplitude (A) and frequency (f) using a control device, as depicted in Figure 4(b). In this study, the vibration was sinusoidal:(7)x=Asin⁡2πft,where x is the displacement of the trough. The real-time vibration was monitored using a laser displacement sensor (Panasonic HG-C1030).In the experiment, a known volume of glass beads (Vtot = 170 mm3) was placed on the middle of the trough surface containing the finlike asperities (Step 1). Then, the trough was vibrated along its central axis, and the glass beads were collected by the silo when the beads were discharged from the discharge hole (+x-direction in Figure 4, Step 2). The collection process was recorded using a digital camera to ascertain the heaping height (hm). The mass of the glass beads was sufficiently small that its influence on the vibration amplitude A was ignored.WhenA=0.3 mm and f = 50 Hz, the heaping height increased over time until all of the glass beads were discharged, as shown in Figure 5. To evaluate the conveyance rate, we selected an interval during which the relationship between hm and t was approximately linear. We considered that the relationship would be linear if the correlation coefficient (r) was greater than 0.9 in this study. The corresponding conveyance rate is(8)Rc,e=πd024Vtot·dhmdt,t∈6,12.Figure 5
Relationship between the experimental heaping height and time.According to the fitting curve in Figure5, Rc,e was 0.089 s−1 in the experiment.In the DEM simulation, an assembly of 1,500 spherical particles (the total occupation was approximatelyVtot = 170 mm3) was employed and completely generated from the particle generation region in 1 s, as shown in Figure 6. The region measured 20 mm × 10 mm. After generation, the trough vibrated along the x-direction for 10 s. The size of the trough and finlike asperity, the inclination angle (θ), and the vibration parameters (A and f) were the same as in the experiment. To shorten the actual time for simulation, we used the particle number instead of the heaping height to define Rc,s:(9)Rc,s=1Ntot·dNcdt,where Ntot = 1,500 is the total particle number and Nc is the discharged particle number. Similarly, Rc,s was also obtained according to the linear part of the t-Nc curve.Figure 6
Sketch of the DEM simulation conditions.The trough and the particle material mechanical parameters are listed in Table1. The particle Poisson’s ratio was based on Tatemoto et al.’s work [38]. The shear modulus was set to 5×107 Pa to improve the computational efficiency. Though this value is small for glass, the influence can be revised using microcontact parameter matching as described below. Other parameters were obtained from the raw material suppliers.Table 1
Material mechanical parameters of the trough and particles.
Parameter Trough Particles Densityρ [kg/m3] 1070 2468 Shear modulusG [Pa] 8.63 × 10 8 5 × 10 7 Poisson’s ratioγ [-] 0.39 0.23In the DEM simulation, the time stepΔt depends primarily on the particle material mechanical parameters [39]. In this study, Δt was 1.51×10-6 s, and the data were recorded every 0.001 s.It is necessary to analyze the significance of the microcontact parameters before matching. To this end, an orthogonal test is simple and appropriate. Generally, the restitution coefficiente and static friction coefficient μs each lies between 0 and 1, and the rolling friction coefficient μr lies between 0 and 0.1. According to Tatemoto et al.’s work [38], the particle-particle restitution coefficient ep-p is 0.9. Assuming that the static and the rolling friction coefficients of the particle-particle and particle-trough interactions are the same (i.e., μs,p-p=μs,p-t=μs and μr,p-p=μr,p-t=μr), the number of the undetermined microcontact parameters is reduced from 5 to 3. The resulting levels of each parameter in the orthogonal test are listed in Table 2.Table 2
Microcontact parameter levels in the orthogonal test.
Level number (A)ep-t [-] (B)μs [-] (C)μr [-] 1 0.3 0.5 0.08 2 0.5 0.3 0.06 3 0.7 0.9 0.04 4 0.9 0.7 0.10The orthogonal test included 3 factors and 4 levels; thus, 16 subsimulations were required. The relationship betweenNc and t in each subsimulation is shown in Figure 7, and the conveyance rates Rc,s are listed in Table 3.Table 3
Orthogonal test scheme and results.
Subsimulation number (A) Level number (B) Level number (C) Level number R c , s (s−1) 1 1 1 1 0.124 2 1 2 2 0.231 3 1 3 3 0.074 4 1 4 4 0.077 5 2 2 4 0.306 6 2 1 3 0.194 7 2 4 2 0.153 8 2 3 1 0.133 9 3 3 2 0.193 10 3 4 1 0.190 11 3 1 4 0.221 12 3 2 3 0.252 13 4 4 3 0.028 14 4 3 4 0.053 15 4 2 1 0.049 16 4 1 2 0.041  RangeRs 0.684 0.390 0.161Figure 7
Relationship between the discharged particle number and vibration time in each subsimulation: (a)ep-t = 0.3, (b) ep-t = 0.5, (c) ep-t = 0.7, and (d) ep-t = 0.9.
(a) (b) (c) (d)The intuitive analysis (range analysis) in the orthogonal tests [40] is a simple method to estimate the significance of the factors. A greater range (Rs) implies that the corresponding factor is more significant. Based on the result in Table 3, we determined that the restitution coefficient of the particle-trough (ep-t) was the most significant factor and that the rolling friction coefficient μr produced less of an effect on the conveyance rate. In addition, the results of simulations numbers 3 and 4 were more similar to Rc,e. Thus, the parameter ep-t was determined to be 0.3 in this study.Although the influence of the rolling friction coefficient on the conveyance rate was not prominent, this influence should be determined as accurately as possible. An effective method to perform this evaluation is the static matching method, also known as AOR matching in this study. The procedure of measuring the AOR is as follows: (a) place a hopper on the horizontal plane; (b) fill the hopper with a certain volume of particles; (c) hoist the hopper slowly until the particles are completely discharged; and (d) measure the AORθr. Figure 8 illustrates the procedure.Figure 8
Procedure to measure the AOR.The experimental value of AOR (θr,e) was 25.4°, while the simulation values (θr,s) based on simulations numbers 3 and 4 were 28.3° and 34.8°, respectively, as shown in Figure 9. The simulated value approached θr,e when μr=0.04 despite μs being larger.Figure 9
AOR: (a)θr,e = 25.4° in the experiment; (b) θr,s = 28.3° in the simulation (number 3); and (c) θr,s = 34.8° in the simulation (number 4).
(a) (b) (c)To determine the final values ofμs and μr, we employed an optimal test whose object was to match (Rc,s,θr,s) with (Rc,e,θr,e), with the conditions(10)Rc,s-Rc,eRc,e≤ξRc∩θr,s-θr,eθr,e≤ξθr,where ξ(Rc) and ξ(θr) were the allowable errors, both set as 0.05. The results are shown in Figure 10, where μr = [0.03, 0.04] and μs = [0.3, 0.5, 0.7, 0.9]. The figure indicates that Rc,s declined as μs rose, while θr,s increased under the same μr. Based on the curves and allowable errors, the optimized intervals of μs (Aop) under various μr were determined. Clearly, Rc was larger when μr = 0.03, and the ultimate values of μs and μr were determined to be 0.8 and 0.03, respectively.Figure 10
Relationship betweenRc,s, θr,s, and μs under various values of μr.By settingep-p=0.9, ep-t=0.3, μs=0.8, and μr=0.03, the conveyance rate Rc,s and the AOR θr,s were 0.091 s−1 and 26.1°, respectively, which further validated the reliability of these predictions. Then, the total microcontact parameters for the DEM simulation were determined. The entire procedure is termed the dynamic-static matching method.


## 3.3. Dimensions of the Finlike Asperity

The parameters to determine the finlike asperity areh, α, and β, as shown in Figure 4. To simplify our exposition, we introduce the following three dimensionless parameters:(11)Yh=hdp,Ya=tan⁡α,Yb=tan⁡β,where dp is the particle diameter (dp=0.5 mm). The levels of each parameter are listed in Table 4. The scope of the present study did not allow simulation of all possible combinations.Table 4
Dimensionless parameters to describe the finlike asperity and associated levels.
Level number Y h Y a Y b 1 0.2 0.1 1 2 0.4 0.2 2 3 0.6 0.3 3 4 0.8 0.4 4 5 1.0 0.5 5Another reason to use the DEM simulation instead of an experiment in this study is that smaller asperities (particularly whenYh < 0.6 and Ya > 0.3) are challenging to manufacture using conventional methods.


## 4. Results and Discussion

### 4.1. Why Do Particles Move Directionally?

Figure11 shows the particle directional conveyance in the trough, where Yh = 0.4, Ya = 0.2, and Yb = 4. A portion of the particles was selected (SP) for analysis. Then, the force acting on the SP was divided into three components: the gravity mpsg, the force in the x-direction Fx, and the force in the z-direction Fz. Fx and Fz were produced by the trough and other particles and satisfied(12)mpsax=Fx-mpsgsin⁡θ,mpsaz=Fz-mpsgcos⁡θ,where mps is the mass of SP and ax and az are the acceleration in the x-direction and z-direction, respectively. Since the directional conveyance occurred in the x-direction, we are concerned only with the first expression in (12).Figure 11
Particle directional conveyance in the trough (Yh = 0.4, Ya = 0.2, and Yb = 4).According to the EDEM postprocessing module,mps was determined to be 1.91×10-5 kg. The relationship between Fx and t is shown in Figure 12(a). If the particles were to vibrate together with the trough, Fx would be symmetric about the axis where Fx=0. However, Fx was asymmetric, and the positive peak force (Fx,u) was greater than the negative peak force (Fx,d) over each vibration period. An asymmetric force was induced by the finlike asperities on the trough surface.Figure 12
(a) Relationship between the force component of the SP along thex-direction and the vibration time. Average velocity components (b) and displacement components (c) of SP and the trough along the x-direction.
(a) (b) (c)The force can alter the object kinestate according to Newton’s laws of motion. The average velocity component of SP along thex-direction (vx) is shown in Figure 12(b). The effect of the asymmetric force Fx, vx was also asymmetric, and its value was generally positive. However, vx of the trough was always symmetric and sinusoidal.Furthermore, the average displacement component of the SP along thex-direction (dx) could be obtained based on vx using(13)dxt=∫0tvxtdt.The relationship betweendx and t is shown in Figure 12(c). The motion of SP could be described as the combination of two submotions. One of the submotions was the local oscillation, and the other was the directional moving toward the +x-direction. The latter submotion was observed to be more effective. The conveyance velocity (vc) was equal to the slope of the envelope line (tan⁡αc). In contrast, the trough oscillated around its equilibrium position only.In fact,vc was not a constant value throughout the conveyance process and was instead related to the packing state of the particles remaining in the trough. This characteristic is generally described by the coordination number (CN), which is defined as the number of particles in contact with a given particle [41]. In other words, if the distance between two particle centers is less than a critical distance, the two particles are considered to be in contact. The critical distance generally ranges from 1.005dp to 1.05dp. This critical distance is suitable for static and quasi-static granular systems in which each particle’s position is nearly unchanged. However, frequent collision occurs among particles in the trough under vibration, implying that particle i has the opportunity to be in contact with particle j at some subsequent time point despite not currently being in contact. Thus, the critical distance can be greater than 1.05dp. In this study, the critical distance is considered as the length of the body diagonal (1.73dp) of a cube whose edge length equals dp, as shown in Figure 13.Figure 13
Critical distance to determine the CN.At any moment, the CN was also not a constant value but ranged from 0 to a maximum value (CNmax). Based on this concept, a dimensionless parameter RCN was introduced to evaluate the packing state and was defined as follows:(14)RCN=NCNNtot-Nc,where NCN was the number of the particles whose CN was greater than np (0≤np≤CNmax). Figure 14 shows the distribution of RCN at t=4 s, at which point the CN ranged from 0 to 10. For a clearer understanding, a median CN (CN50) was introduced, defined as the CN at which the corresponding RCN equals 0.5 in the CN distribution (np-RCN) curve. A larger CN50 indicates that the particles are packed together more closely.Figure 14
(a) The CN distribution of particles in the trough and (b) the CN distribution curve and the median CN (CN50).
(a) (b)Throughout the simulation, the relationship betweendx and t was tracked and is shown in Figure 15. As the vibration time t increased, the particles were gradually discharged from the trough, and CN50 subsequently declined. In contrast, vc increased until all of the particles were discharged. However, vc did not change sharply with t but rather maintained a steady value over an extended time. For example, vc=11.85 mm/s from t=1.5 s to 6.5 s under this condition. Therefore, instead of Rc, this steady vc (vcs) served as a more intuitive parameter to evaluate the particle directional conveyance performance in DEM simulations.Figure 15
Particles average displacement in thex-direction throughout the simulation and the variation of the CN distribution and vc at various time points (t=0.5 s, t=4 s, and t=7 s).v c at time point t was obtained by calculating the slope of the envelope line containing three adjacent vibration periods (T) in the t-dx curve, as shown in Figure 12(c).


### 4.2. Influence of the Asperity Height on the Conveyance Process

WithYa=0.2 and Yb=4, the resulting relationship between vcs and Yh is shown in Figure 16. When Yh=0.2 and 0.4, vcs was nearly the same. When Yh was greater than 0.4, vcs declined as Yh increased. Under this condition, the finlike asperity shape was the same but with proportional size, implying that a greater Yh results in a larger asperity and that the distance between two asperities (da) also increases. In this work, we introduce the parameter φa to quantitatively describe the relationship between da and dp, as the ratio of da to dp:(15)φa=dadp=Yh1Ya+1Yb.Since the vibration parameters were invariable, the energy provided for the particles was constant and finite. Thus, it was less likely for the particles to leap over the asperities with larger da. According to the results in Figure 16, the directional conveyance was weakened when φa was greater than approximately 2.Figure 16
The influence of the finlike asperity height on the particle conveyance velocity and the distance between two asperities (Ya = 0.2 and Yb = 4).


### 4.3. Influence of the Asperity Slope Inclination Angles on the Conveyance Process

KeepingYh=0.4 and Yb=4, the resulting relationships between Ya and vcs and Ya and φa are shown in Figure 17(a). The figure indicates that as the asperity slope inclination angle α increased, the distance between two asperities da decreased monotonously. However, the conveyance velocity vcs varied circuitously. When the coefficient Ya was less than 0.3, that is, when φa was greater than 1.43, vcs increased with increasing Ya. When Ya was greater than 0.3, vcs instead was reduced, and the particles even flowed toward the opposite direction when Ya = 0.5. This result also shows that a much smaller da does not contribute to the particle directional conveyance. A favorable range of φa under this condition was from 1.1 to 2.1.Figure 17
Influence of the asperity slope inclination angle on the conveyance velocity and the distance between two asperities. (a) Influence ofα when Yh=0.4 and Yb=4 and (b) influence of β when Yh=0.4 and Ya=0.2.
(a) (b)KeepingYh=0.4 and Ya=0.2, the resulting relationships between Yb and vcs and Yb and φa are shown in Figure 17(b). The figure indicates that da declined with increasing β. However, since β > 45°, Yb > 1, and the effect on da was small relative to the effect on α. Therefore, vcs did not change substantially either, particularly when Yb≥2.A familiar configuration of the finlike asperity is whenα+β=90°, that is, when the finlike asperity shape is a right triangle. By defining Yb=1/Ya and Ya = [0.2, 0.25, 0.3, 0.35, 0.4], we ensure that Yb>2. As β was not a sensitive parameter to vcs when Yb≥2, the focus of this work was the effects of h and α. In addition, Yh was set as 0.4, 0.5, 0.6, 0.7, and 0.8 for a denser micromesh analysis. The relationships between Ya and vcs and Ya and φa under various Yh are shown in Figure 18. The curves in Figure 18(a) indicate that with increasing h, the influence of α on vcs weakened gradually, and h had more of a significant influence than α within this scope.Figure 18
(a) The influence ofYh and Ya on vcs and (b) the effects of Yh and Ya on φa.
(a) (b)Figure18(b) reconfirms that da declined with increasing α and increased as h rose. Correspondingly, when φa ranged from 1 to 2, vcs initially increased and then decreased. When φa decreased to 2, vcs increased monotonously, indicating that the particle conveyance performance was improved within 1<φa<2. However, da could not be used to accurately evaluate vcs. For example, when Yh=0.4 and Ya=0.25 and when Yh=0.6 and Ya=0.4, φa was almost 1.7, but vcs was 14.75 mm/s and 7.15 mm/s, respectively. Despite the fact that da is not suitable to quantitatively evaluate vcs, this parameter can provide guidance to roughly optimize h and α to improve the particle conveyance performance.


### 4.4. Discussion

According to the DEM simulation, the particle directional conveyance under vibration is due to the asymmetric forces in the longitudinal direction. This mechanism is similar to material conveyance in a nonharmonic vibratory conveyor [10]. However, the asymmetric force is caused by the finlike asperities in our study rather than the trough’s nonsinusoidal vibration in the traditional vibratory conveyance. Using the trough with finlike asperities, the vibration form can be sinusoidal, which benefits the simplification of the driving mechanism. This simplification shows great significance for nonterrestrial tasks.CN is an appropriate parameter to evaluate the influence of the particle packing state on the conveyance capacity. In particular, the CN characterizes the advantage over the particle layer thickness [42] when there is only one layer. Of course, the use of CN is feasible only in simulation studies and not experimental studies. The distribution of CN at a given time point is described by the median CN, CN50. A large CN50 indicates that more particles are or will be in contact, which can reduce directional conveyance.The simulation results also show the influence of the finlike asperity heighth and slope inclination angles α and β on the particle directional conveyance capacity. The values of h, α, and β determine the asperity interval da. An excessively small or large da does not benefit particle directional conveyance, and an intermediately sized da (1dp~2dp) is more appropriate. The importance of the right asperity size is also validated in Farkas et al.’ study despite the vibration direction being vertical [25], which indicates that the finlike asperity size should match the particle diameter dp to improve the conveyance capacity.


## 4.1. Why Do Particles Move Directionally?

Figure11 shows the particle directional conveyance in the trough, where Yh = 0.4, Ya = 0.2, and Yb = 4. A portion of the particles was selected (SP) for analysis. Then, the force acting on the SP was divided into three components: the gravity mpsg, the force in the x-direction Fx, and the force in the z-direction Fz. Fx and Fz were produced by the trough and other particles and satisfied(12)mpsax=Fx-mpsgsin⁡θ,mpsaz=Fz-mpsgcos⁡θ,where mps is the mass of SP and ax and az are the acceleration in the x-direction and z-direction, respectively. Since the directional conveyance occurred in the x-direction, we are concerned only with the first expression in (12).Figure 11
Particle directional conveyance in the trough (Yh = 0.4, Ya = 0.2, and Yb = 4).According to the EDEM postprocessing module,mps was determined to be 1.91×10-5 kg. The relationship between Fx and t is shown in Figure 12(a). If the particles were to vibrate together with the trough, Fx would be symmetric about the axis where Fx=0. However, Fx was asymmetric, and the positive peak force (Fx,u) was greater than the negative peak force (Fx,d) over each vibration period. An asymmetric force was induced by the finlike asperities on the trough surface.Figure 12
(a) Relationship between the force component of the SP along thex-direction and the vibration time. Average velocity components (b) and displacement components (c) of SP and the trough along the x-direction.
(a) (b) (c)The force can alter the object kinestate according to Newton’s laws of motion. The average velocity component of SP along thex-direction (vx) is shown in Figure 12(b). The effect of the asymmetric force Fx, vx was also asymmetric, and its value was generally positive. However, vx of the trough was always symmetric and sinusoidal.Furthermore, the average displacement component of the SP along thex-direction (dx) could be obtained based on vx using(13)dxt=∫0tvxtdt.The relationship betweendx and t is shown in Figure 12(c). The motion of SP could be described as the combination of two submotions. One of the submotions was the local oscillation, and the other was the directional moving toward the +x-direction. The latter submotion was observed to be more effective. The conveyance velocity (vc) was equal to the slope of the envelope line (tan⁡αc). In contrast, the trough oscillated around its equilibrium position only.In fact,vc was not a constant value throughout the conveyance process and was instead related to the packing state of the particles remaining in the trough. This characteristic is generally described by the coordination number (CN), which is defined as the number of particles in contact with a given particle [41]. In other words, if the distance between two particle centers is less than a critical distance, the two particles are considered to be in contact. The critical distance generally ranges from 1.005dp to 1.05dp. This critical distance is suitable for static and quasi-static granular systems in which each particle’s position is nearly unchanged. However, frequent collision occurs among particles in the trough under vibration, implying that particle i has the opportunity to be in contact with particle j at some subsequent time point despite not currently being in contact. Thus, the critical distance can be greater than 1.05dp. In this study, the critical distance is considered as the length of the body diagonal (1.73dp) of a cube whose edge length equals dp, as shown in Figure 13.Figure 13
Critical distance to determine the CN.At any moment, the CN was also not a constant value but ranged from 0 to a maximum value (CNmax). Based on this concept, a dimensionless parameter RCN was introduced to evaluate the packing state and was defined as follows:(14)RCN=NCNNtot-Nc,where NCN was the number of the particles whose CN was greater than np (0≤np≤CNmax). Figure 14 shows the distribution of RCN at t=4 s, at which point the CN ranged from 0 to 10. For a clearer understanding, a median CN (CN50) was introduced, defined as the CN at which the corresponding RCN equals 0.5 in the CN distribution (np-RCN) curve. A larger CN50 indicates that the particles are packed together more closely.Figure 14
(a) The CN distribution of particles in the trough and (b) the CN distribution curve and the median CN (CN50).
(a) (b)Throughout the simulation, the relationship betweendx and t was tracked and is shown in Figure 15. As the vibration time t increased, the particles were gradually discharged from the trough, and CN50 subsequently declined. In contrast, vc increased until all of the particles were discharged. However, vc did not change sharply with t but rather maintained a steady value over an extended time. For example, vc=11.85 mm/s from t=1.5 s to 6.5 s under this condition. Therefore, instead of Rc, this steady vc (vcs) served as a more intuitive parameter to evaluate the particle directional conveyance performance in DEM simulations.Figure 15
Particles average displacement in thex-direction throughout the simulation and the variation of the CN distribution and vc at various time points (t=0.5 s, t=4 s, and t=7 s).v c at time point t was obtained by calculating the slope of the envelope line containing three adjacent vibration periods (T) in the t-dx curve, as shown in Figure 12(c).


## 4.2. Influence of the Asperity Height on the Conveyance Process

WithYa=0.2 and Yb=4, the resulting relationship between vcs and Yh is shown in Figure 16. When Yh=0.2 and 0.4, vcs was nearly the same. When Yh was greater than 0.4, vcs declined as Yh increased. Under this condition, the finlike asperity shape was the same but with proportional size, implying that a greater Yh results in a larger asperity and that the distance between two asperities (da) also increases. In this work, we introduce the parameter φa to quantitatively describe the relationship between da and dp, as the ratio of da to dp:(15)φa=dadp=Yh1Ya+1Yb.Since the vibration parameters were invariable, the energy provided for the particles was constant and finite. Thus, it was less likely for the particles to leap over the asperities with larger da. According to the results in Figure 16, the directional conveyance was weakened when φa was greater than approximately 2.Figure 16
The influence of the finlike asperity height on the particle conveyance velocity and the distance between two asperities (Ya = 0.2 and Yb = 4).


## 4.3. Influence of the Asperity Slope Inclination Angles on the Conveyance Process

KeepingYh=0.4 and Yb=4, the resulting relationships between Ya and vcs and Ya and φa are shown in Figure 17(a). The figure indicates that as the asperity slope inclination angle α increased, the distance between two asperities da decreased monotonously. However, the conveyance velocity vcs varied circuitously. When the coefficient Ya was less than 0.3, that is, when φa was greater than 1.43, vcs increased with increasing Ya. When Ya was greater than 0.3, vcs instead was reduced, and the particles even flowed toward the opposite direction when Ya = 0.5. This result also shows that a much smaller da does not contribute to the particle directional conveyance. A favorable range of φa under this condition was from 1.1 to 2.1.Figure 17
Influence of the asperity slope inclination angle on the conveyance velocity and the distance between two asperities. (a) Influence ofα when Yh=0.4 and Yb=4 and (b) influence of β when Yh=0.4 and Ya=0.2.
(a) (b)KeepingYh=0.4 and Ya=0.2, the resulting relationships between Yb and vcs and Yb and φa are shown in Figure 17(b). The figure indicates that da declined with increasing β. However, since β > 45°, Yb > 1, and the effect on da was small relative to the effect on α. Therefore, vcs did not change substantially either, particularly when Yb≥2.A familiar configuration of the finlike asperity is whenα+β=90°, that is, when the finlike asperity shape is a right triangle. By defining Yb=1/Ya and Ya = [0.2, 0.25, 0.3, 0.35, 0.4], we ensure that Yb>2. As β was not a sensitive parameter to vcs when Yb≥2, the focus of this work was the effects of h and α. In addition, Yh was set as 0.4, 0.5, 0.6, 0.7, and 0.8 for a denser micromesh analysis. The relationships between Ya and vcs and Ya and φa under various Yh are shown in Figure 18. The curves in Figure 18(a) indicate that with increasing h, the influence of α on vcs weakened gradually, and h had more of a significant influence than α within this scope.Figure 18
(a) The influence ofYh and Ya on vcs and (b) the effects of Yh and Ya on φa.
(a) (b)Figure18(b) reconfirms that da declined with increasing α and increased as h rose. Correspondingly, when φa ranged from 1 to 2, vcs initially increased and then decreased. When φa decreased to 2, vcs increased monotonously, indicating that the particle conveyance performance was improved within 1<φa<2. However, da could not be used to accurately evaluate vcs. For example, when Yh=0.4 and Ya=0.25 and when Yh=0.6 and Ya=0.4, φa was almost 1.7, but vcs was 14.75 mm/s and 7.15 mm/s, respectively. Despite the fact that da is not suitable to quantitatively evaluate vcs, this parameter can provide guidance to roughly optimize h and α to improve the particle conveyance performance.


## 4.4. Discussion

According to the DEM simulation, the particle directional conveyance under vibration is due to the asymmetric forces in the longitudinal direction. This mechanism is similar to material conveyance in a nonharmonic vibratory conveyor [10]. However, the asymmetric force is caused by the finlike asperities in our study rather than the trough’s nonsinusoidal vibration in the traditional vibratory conveyance. Using the trough with finlike asperities, the vibration form can be sinusoidal, which benefits the simplification of the driving mechanism. This simplification shows great significance for nonterrestrial tasks.CN is an appropriate parameter to evaluate the influence of the particle packing state on the conveyance capacity. In particular, the CN characterizes the advantage over the particle layer thickness [42] when there is only one layer. Of course, the use of CN is feasible only in simulation studies and not experimental studies. The distribution of CN at a given time point is described by the median CN, CN50. A large CN50 indicates that more particles are or will be in contact, which can reduce directional conveyance.The simulation results also show the influence of the finlike asperity heighth and slope inclination angles α and β on the particle directional conveyance capacity. The values of h, α, and β determine the asperity interval da. An excessively small or large da does not benefit particle directional conveyance, and an intermediately sized da (1dp~2dp) is more appropriate. The importance of the right asperity size is also validated in Farkas et al.’ study despite the vibration direction being vertical [25], which indicates that the finlike asperity size should match the particle diameter dp to improve the conveyance capacity.


## 5. Conclusions

For the purpose of granular material transport on future lunar base construction, a type of vibratory conveyance method that considers finlike asperities on a trough surface was introduced. In this study, we employed the DEM to investigate the particle directional conveyance characteristics in the trough. The following conclusions were drawn.(a)
The restitution coefficient between the trough and the particleep-t is the factor that most strongly influences the particle directional conveyance capacity. The two most important factors are the static friction coefficient μs and the rolling friction coefficient μr.(b)
Due to the finlike asperities, particles are exposed to an asymmetric forceFx under longitudinal and sinusoidal vibrations. The cumulative effect of Fx over time governs the particle directional conveyance. During the conveyance process, the median CN, CN50, affects the conveyance velocity vc; a greater CN50 induces a lower vc because the multitude of contacting particle collisions weakens the effects of conveyance.(c)
The particle steady conveyance velocityvcs is affected by the asperity height h and the slope inclination angles α and β. In fact, these parameters directly determine the distance between the two finlike asperities da, and an undersized or oversized da tends to reduce vcs. In this study, vcs increased when da was between one and two times the particle diameter dp.



---

*Source: 8260462-2018-02-26.xml*