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.
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.
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:where and are the translational and angular velocity vectors of particle , respectively; and are the mass and the rotational inertia of particle , respectively; is the gravitational acceleration; is the number of particles in collision with particle ; and and are the contact force and torque generated between particles and , respectively, and are composed of several components: namely,where and are the normal and the tangential forces, respectively, and and are the torques caused by the normal force and the tangential force, respectively. 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 and . 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 Tangential force Torque
The parameters and are the static and rolling friction coefficients, respectively; and are the normal overlap and tangential displacements, respectively; and are the relative normal and tangential velocities, respectively; is the unit vector pointing from the center of particle to particle ; is the unit angular velocity vector when a relative rotation occurs; and , , and are the equivalent Young’s modulus, particle radius, and mass, respectively [32].
In (4a), (4b), and (4c), the superscript is the time step number, is the elastic component of the tangential force, and and are the tangential displacement increment and the normal force increment, respectively. depends on the interaction phase. If the interaction is in the loading phase, the increment of is simply proportional to (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 parameters , , and are given bywhere is the equivalent shear modulus and 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 , , and , 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.
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 (), 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 4 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 (), as shown in Figure 4(a). The trough was positioned obliquely to obtain an inclination angle () of 10°. The height of each asperity () 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 () 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.
(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 () and frequency () using a control device, as depicted in Figure 4(b). In this study, the vibration was sinusoidal:where 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 ( = 170 mm3) was placed on the middle of the trough surface containing the finlike asperities (Step ). 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 (+-direction in Figure 4, Step ). The collection process was recorded using a digital camera to ascertain the heaping height (). The mass of the glass beads was sufficiently small that its influence on the vibration amplitude was ignored.
When mm and = 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 and was approximately linear. We considered that the relationship would be linear if the correlation coefficient () was greater than 0.9 in this study. The corresponding conveyance rate is
According to the fitting curve in Figure 5, was 0.089 s−1 in the experiment.
In the DEM simulation, an assembly of 1,500 spherical particles (the total occupation was approximately = 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 -direction for 10 s. The size of the trough and finlike asperity, the inclination angle (), and the vibration parameters ( and ) 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 :where = 1,500 is the total particle number and is the discharged particle number. Similarly, was also obtained according to the linear part of the - curve.
The trough and the particle material mechanical parameters are listed in Table 1. The particle Poisson’s ratio was based on Tatemoto et al.’s work [38]. The shear modulus was set to 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.
In the DEM simulation, the time step depends primarily on the particle material mechanical parameters [39]. In this study, was 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 coefficient and static friction coefficient each lies between 0 and 1, and the rolling friction coefficient lies between 0 and 0.1. According to Tatemoto et al.’s work [38], the particle-particle restitution coefficient 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., and ), 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.
The orthogonal test included 3 factors and 4 levels; thus, 16 subsimulations were required. The relationship between and in each subsimulation is shown in Figure 7, and the conveyance rates are listed in Table 3.
(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 () 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 () was the most significant factor and that the rolling friction coefficient produced less of an effect on the conveyance rate. In addition, the results of simulations numbers 3 and 4 were more similar to . Thus, the parameter 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 . Figure 8 illustrates the procedure.
The experimental value of AOR () was 25.4°, while the simulation values () based on simulations numbers 3 and 4 were 28.3° and 34.8°, respectively, as shown in Figure 9. The simulated value approached when despite being larger.
(a)
(b)
(c)
To determine the final values of and , we employed an optimal test whose object was to match () with (), with the conditionswhere and were the allowable errors, both set as 0.05. The results are shown in Figure 10, where = [0.03, 0.04] and = [0.3, 0.5, 0.7, 0.9]. The figure indicates that declined as rose, while increased under the same . Based on the curves and allowable errors, the optimized intervals of () under various were determined. Clearly, was larger when = 0.03, and the ultimate values of and were determined to be 0.8 and 0.03, respectively.
By setting , , , and , the conveyance rate and the AOR 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 are , , and , as shown in Figure 4. To simplify our exposition, we introduce the following three dimensionless parameters:where is the particle diameter ( 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.
Another reason to use the DEM simulation instead of an experiment in this study is that smaller asperities (particularly when < 0.6 and > 0.3) are challenging to manufacture using conventional methods.
4. Results and Discussion
4.1. Why Do Particles Move Directionally?
Figure 11 shows the particle directional conveyance in the trough, where = 0.4, = 0.2, and = 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 , the force in the -direction , and the force in the -direction . and were produced by the trough and other particles and satisfiedwhere is the mass of SP and and are the acceleration in the -direction and -direction, respectively. Since the directional conveyance occurred in the -direction, we are concerned only with the first expression in (12).
According to the EDEM postprocessing module, was determined to be kg. The relationship between and is shown in Figure 12(a). If the particles were to vibrate together with the trough, would be symmetric about the axis where . However, was asymmetric, and the positive peak force () was greater than the negative peak force () over each vibration period. An asymmetric force was induced by the finlike asperities on the trough surface.
(a)
(b)
(c)
The force can alter the object kinestate according to Newton’s laws of motion. The average velocity component of SP along the -direction () is shown in Figure 12(b). The effect of the asymmetric force , was also asymmetric, and its value was generally positive. However, of the trough was always symmetric and sinusoidal.
Furthermore, the average displacement component of the SP along the -direction () could be obtained based on using
The relationship between and 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 +-direction. The latter submotion was observed to be more effective. The conveyance velocity () was equal to the slope of the envelope line (). In contrast, the trough oscillated around its equilibrium position only.
In fact, 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.005 to 1.05. 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 has the opportunity to be in contact with particle at some subsequent time point despite not currently being in contact. Thus, the critical distance can be greater than 1.05. In this study, the critical distance is considered as the length of the body diagonal (1.73) of a cube whose edge length equals , as shown in Figure 13.
At any moment, the CN was also not a constant value but ranged from 0 to a maximum value (). Based on this concept, a dimensionless parameter was introduced to evaluate the packing state and was defined as follows:where was the number of the particles whose CN was greater than (). Figure 14 shows the distribution of at 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 equals 0.5 in the CN distribution () curve. A larger CN50 indicates that the particles are packed together more closely.
(a)
(b)
Throughout the simulation, the relationship between and was tracked and is shown in Figure 15. As the vibration time increased, the particles were gradually discharged from the trough, and CN50 subsequently declined. In contrast, increased until all of the particles were discharged. However, did not change sharply with but rather maintained a steady value over an extended time. For example, mm/s from s to 6.5 s under this condition. Therefore, instead of , this steady () served as a more intuitive parameter to evaluate the particle directional conveyance performance in DEM simulations.
at time point was obtained by calculating the slope of the envelope line containing three adjacent vibration periods () in the - curve, as shown in Figure 12(c).
4.2. Influence of the Asperity Height on the Conveyance Process
With and , the resulting relationship between and is shown in Figure 16. When and 0.4, was nearly the same. When was greater than 0.4, declined as increased. Under this condition, the finlike asperity shape was the same but with proportional size, implying that a greater results in a larger asperity and that the distance between two asperities () also increases. In this work, we introduce the parameter to quantitatively describe the relationship between and , as the ratio of to :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 . According to the results in Figure 16, the directional conveyance was weakened when was greater than approximately 2.
4.3. Influence of the Asperity Slope Inclination Angles on the Conveyance Process
Keeping and , the resulting relationships between and and and are shown in Figure 17(a). The figure indicates that as the asperity slope inclination angle increased, the distance between two asperities decreased monotonously. However, the conveyance velocity varied circuitously. When the coefficient was less than 0.3, that is, when was greater than 1.43, increased with increasing . When was greater than 0.3, instead was reduced, and the particles even flowed toward the opposite direction when = 0.5. This result also shows that a much smaller does not contribute to the particle directional conveyance. A favorable range of under this condition was from 1.1 to 2.1.
(a)
(b)
Keeping and , the resulting relationships between and and and are shown in Figure 17(b). The figure indicates that declined with increasing . However, since > 45°, > 1, and the effect on was small relative to the effect on . Therefore, did not change substantially either, particularly when .
A familiar configuration of the finlike asperity is when , that is, when the finlike asperity shape is a right triangle. By defining and = [0.2, 0.25, 0.3, 0.35, 0.4], we ensure that . As was not a sensitive parameter to when , the focus of this work was the effects of and . In addition, was set as 0.4, 0.5, 0.6, 0.7, and 0.8 for a denser micromesh analysis. The relationships between and and and under various are shown in Figure 18. The curves in Figure 18(a) indicate that with increasing , the influence of on weakened gradually, and had more of a significant influence than within this scope.
(a)
(b)
Figure 18(b) reconfirms that declined with increasing and increased as rose. Correspondingly, when ranged from 1 to 2, initially increased and then decreased. When decreased to 2, increased monotonously, indicating that the particle conveyance performance was improved within . However, could not be used to accurately evaluate . For example, when and and when and , was almost 1.7, but was 14.75 mm/s and 7.15 mm/s, respectively. Despite the fact that is not suitable to quantitatively evaluate , this parameter can provide guidance to roughly optimize 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 height and slope inclination angles and on the particle directional conveyance capacity. The values of , , and determine the asperity interval . An excessively small or large does not benefit particle directional conveyance, and an intermediately sized (12) 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 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 particle is the factor that most strongly influences the particle directional conveyance capacity. The two most important factors are the static friction coefficient and the rolling friction coefficient .(b)Due to the finlike asperities, particles are exposed to an asymmetric force under longitudinal and sinusoidal vibrations. The cumulative effect of over time governs the particle directional conveyance. During the conveyance process, the median CN, CN50, affects the conveyance velocity ; a greater CN50 induces a lower because the multitude of contacting particle collisions weakens the effects of conveyance.(c)The particle steady conveyance velocity is affected by the asperity height and the slope inclination angles and . In fact, these parameters directly determine the distance between the two finlike asperities , and an undersized or oversized tends to reduce . In this study, increased when was between one and two times the particle diameter .
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
The authors are grateful for the support of the National Natural Science Foundation of China (Grant no. 51575122) and the Self-Planned Task of State Key Laboratory of Robotics and System (HIT) (Grant no. SKLRS201616B).