DEM Study on the Segregation of a Non-Spherical Intruder in a Vibrated Granular Bed

: The segregation process of a single large intruder in a vibrated bed of small particles has been widely studied, but most previous studies focused on spherical intruders. In this work, the discrete element method was used to study the effects of vibration conditions and intruder shape on the dimensionless ascending velocity ( v a ) of the intruder. The intruder was in a prolate shape with aspect ratio varied but its equivalent diameter ﬁxed. Three equivalent diameters, namely volume-equivalent diameter, surface-area-equivalent diameter, and Sauter diameter, were used. It was found that v a increases and then decreases with the rise of the dimensionless vibration amplitude ( A d ) and the dimensionless vibration frequency ( f d ), and v a increases with the decrease of the sphericity of the intruder ( Φ ). Moreover, the porosity variation in the vibrated bed and the granular temperature were analyzed, which can be linked to the change of v a . It was further found that v a can be uniformly correlated to A d · f 0.5 d , while the critical change of the response of v a to A d and f d occurs at Γ = 4.83, where Γ is the vibration intensity. Based on these ﬁndings, a piecewise equation was proposed to predict v a as a function of A d , f d , and Φ .


Introduction
The segregation of granular mixtures under vibration is often encountered in various industrial processes [1][2][3]. A good understanding of the segregation mechanism can help the optimization and control of the related processes. The research of such a phenomenon often starts with the segregation of a single large intruder in an otherwise homogeneous granular bed of small particles [4][5][6]. Under vibration, the large intruder normally ascends in the granular bed. Based on the statistical analysis on the interaction between the intruder and the surrounding small particles, several kinds of segregation mechanisms have been proposed for such an ascending phenomenon, including void filling, global convection, etc. [7][8][9][10][11]. In addition, in different experimental and numerical studies, the ascending velocity, which is the average velocity of the intruder in its rising process from the bottom to the top of the vibrated bed, was often modeled against different controlling variables [12].
The literature indicates the dependence of the ascending velocity on the following variables: particle properties including the density ratio (ρ r ) [13][14][15] and the size ratio (d r ) [16,17]; particle bed features including the aspect ratio [18], friction [19], and filling height [20]; and vibration conditions including vibration amplitude [21], frequency [12], and intensity [10,22]. For particle properties, a general observation is that the increase of the size ratio and density ratio both increase the ascending velocity [17,20,22]. The shape of the container would affect the granular flow direction [18,23] and control the particle convection [24]. In addition, with or without friction between particles and the wall, the rise time of the intruder is rather different [25]. For vibration conditions, the increase of Table 1. List of models for the ascending of a spherical intruder in a vibrated bed of small particles.
In this paper, the segregation process of a single non-spherical intruder is studied by the discrete element method (DEM). The effects of vibration amplitude and frequency on the defined dimensionless ascending velocity are investigated. Based on the simulation results, an equation is proposed to uniformly link the ascending velocity to the vibration conditions and the shape of the intruder, and particularly the non-monotonic change of the ascending velocity with the vibration conditions is also uniformly considered. These results provide new insights into particle segregation under vibration and new dimensionless models to predict the segregation speed, which could be useful for modeling powder mixing and segregation.

Model Description
The discrete element method (DEM) has been proven as an effective technology to study granular mixtures at a particle scale [34][35][36][37]. Considering computational efficiency and successful experience according to previous literature, the typical Hertz-Mindlin particle contact model is adopted in this work. Based on Newton s laws of motion, the translational motion and rotational motion of particle i are, respectively, determined by Processes 2021, 9, 448 3 of 15 where v i , ω i , and I i are the translational and angular velocities, and the moment of inertia of particle i, respectively; g is the gravitational acceleration and t is time; m i is mass of particle I; n i is the total number of particles in contact with particle i; F n,ij is the normal contact force, including the normal-contact force and the normal-damping force. F t,ij is the tangential contact force, including the tangential contact and the tangential damping force. The commercial software EDEM 2018 (DEM Solutions Ltd., Edinburgh, United Kingdom) was used to conduct the DEM simulations, and the details of the force models can be found in the manual [38].

Simulation Conditions and Validation
A 3D cuboid "container" without physical walls was adopted in this research, as shown in Figure 1a. The bed size was 20 × 20 × 60 mm and periodic boundary conditions (PBCs) were applied along the X-axis and Y-axis directions, which can prevent granular convection induced by the side wall during the vibration [24,39] and its effect on the ascending of the intruder [29]. The initial undisturbed bed consisted of one single nonspherical large intruder in an otherwise homogeneous spherical granular bed with the total bed height H bed ≈ 30·d s , where d s is the diameter of small particles. The parameters used in the simulations are given in Table 2. Continuous and sinusoidal vibrations with wide range of frequency f and amplitude A were applied to the bottom wall along the Z-axis. The non-spherical intruder was in a prolate shape and was initially placed with its long axis aligning with the Z-axis direction. According to our observations, the orientation of the intruder had slight changes during the ascending process, and the average angle between its long axis and the vibration direction was less than 15 degrees, as shown in Figure 1b.
It should be noted that there are different ways to consider the size of a non-spherical particle. In the present research, the equivalent diameter of a non-spherical intruder was calculated in three ways; namely, volume equivalent diameter (d Vol ), area equivalent diameter (d Area ), and Sauter diameter (d Sau ), which are respectively given by where V and S are, respectively, the volume and surface of the non-spherical intruder. The axial sizes of the non-spherical intruders of different sphericities are listed in Table 3.
where vi, ωi, and Ii are the translational and angular velocities, and the moment of inertia of particle i, respectively; g is the gravitational acceleration and t is time; mi is mass of particle I; ni is the total number of particles in contact with particle i; Fn,ij is the normal contact force, including the normal-contact force and the normal-damping force. Ft,ij is the tangential contact force, including the tangential contact and the tangential damping force. The commercial software EDEM 2018 (DEM Solutions Ltd., Edinburgh, United Kingdom) was used to conduct the DEM simulations, and the details of the force models can be found in the manual [38].

Simulation Conditions and Validation
A 3D cuboid "container" without physical walls was adopted in this research, as shown in Figure 1a. The bed size was 20 × 20 × 60 mm and periodic boundary conditions (PBCs) were applied along the X-axis and Y-axis directions, which can prevent granular convection induced by the side wall during the vibration [24,39] and its effect on the ascending of the intruder [29]. The initial undisturbed bed consisted of one single nonspherical large intruder in an otherwise homogeneous spherical granular bed with the total bed height Hbed ≈ 30·ds, where ds is the diameter of small particles. The parameters used in the simulations are given in Table 2. Continuous and sinusoidal vibrations with wide range of frequency f and amplitude A were applied to the bottom wall along the Zaxis. The non-spherical intruder was in a prolate shape and was initially placed with its long axis aligning with the Z-axis direction. According to our observations, the orientation of the intruder had slight changes during the ascending process, and the average angle between its long axis and the vibration direction was less than 15 degrees, as shown in Figure 1b.

Value of Parameter
Poisson's ratio of particles 0.25 Poisson's ratio of container base 0.29 Young's modulus of particles (Pa) 2.01 × 10 5 Young's modulus of container base (Pa) 2.55 × 10 10 Coefficient of restitution: particle-particle 0.5 Coefficient of static friction: particle-particle 0.5 Coefficient of rolling friction: particle-particle 0.01 Coefficient of restitution: particle-base 0.5 Coefficient of static friction: particle-base 0.5 Coefficient of rolling friction: particle-base 0.01 Intruder equivalent diameter, d l (mm) 6  The simulation time step was set to 9.06 × 10 −5 s, which is 20% of the Rayleigh time step. The simulation time for each case was over 25 s. To validate the model, the simulated results were compared to the literature. The simulated ascending process is shown in Figure 2. It can be seen that the height-time curve simulated by our model is in good agreement with that of the physical experiments [11], as shown in Figure 2a. In addition, as shown in Figure 2b, the simulated total ascending time as a function of the vibration amplitude is in good agreement with the experiment results [11]. The good agreement can validate the DEM model. results were compared to the literature. The simulated ascending process is shown in Figure 2. It can be seen that the height-time curve simulated by our model is in good agreement with that of the physical experiments [11], as shown in Figure 2a. In addition, as shown in Figure 2b, the simulated total ascending time as a function of the vibration amplitude is in good agreement with the experiment results [11]. The good agreement can validate the DEM model.

Evaluation Indexes
In this research, the ascending velocity of the intruder, va, is defined as the average vertical displacement of the intruder per oscillation cycle [29]: where ΔH is the vertical distance traveled by the intruder in the ascending process, mm; ds is the size of small particles; T is the number of oscillation cycles in the ascending process. It is worth noting that va obtained in Equation (6) is an average value but is reasonable to represent the ascending behavior. As shown in Figure 3, with different vibration amplitudes, the ascending height is generally in linear dependence to the vibration time, and the fluctuations are observed only when the height of the intruder is close to the top of the bed.

Evaluation Indexes
In this research, the ascending velocity of the intruder, v a , is defined as the average vertical displacement of the intruder per oscillation cycle [29]: where ∆H is the vertical distance traveled by the intruder in the ascending process, mm; d s is the size of small particles; T is the number of oscillation cycles in the ascending process.
It is worth noting that v a obtained in Equation (6) is an average value but is reasonable to represent the ascending behavior. As shown in Figure 3, with different vibration amplitudes, the ascending height is generally in linear dependence to the vibration time, and the fluctuations are observed only when the height of the intruder is close to the top of the bed. ure 2. It can be seen that the height-time curve simulated by our model is in good agreement with that of the physical experiments [11], as shown in Figure 2a. In addition, as shown in Figure 2b, the simulated total ascending time as a function of the vibration amplitude is in good agreement with the experiment results [11]. The good agreement can validate the DEM model.

Evaluation Indexes
In this research, the ascending velocity of the intruder, va, is defined as the average vertical displacement of the intruder per oscillation cycle [29]: where ΔH is the vertical distance traveled by the intruder in the ascending process, mm; ds is the size of small particles; T is the number of oscillation cycles in the ascending process. It is worth noting that va obtained in Equation (6) is an average value but is reasonable to represent the ascending behavior. As shown in Figure 3, with different vibration amplitudes, the ascending height is generally in linear dependence to the vibration time, and the fluctuations are observed only when the height of the intruder is close to the top of the bed.

Effects of Vibration Amplitude on v a
As shown in Figure 4a,b, v a of a non-spherical intruder increased first and then decreased with dimensionless vibration amplitude A d (A d = A/d s ), and the critical change occurred at A d = 1.5. In addition, v a of the non-spherical intruder of aspect ratio a/b = 3.0 was larger than that of a/b = 2.0. Additionally, regardless of which equivalent diameter was used, the non-spherical intruder's ascending velocity was much higher than that of a

Effects of Vibration Amplitude on va
As shown in Figure 4a,b, va of a non-spherical intruder increased first and then decreased with dimensionless vibration amplitude Ad (Ad = A/ds), and the critical change occurred at Ad = 1.5. In addition, va of the non-spherical intruder of aspect ratio a/b = 3.0 was larger than that of a/b = 2.0. Additionally, regardless of which equivalent diameter was used, the non-spherical intruder's ascending velocity was much higher than that of a spherical intruder (a/b = 1.0), and the velocity difference increased with the increase of the aspect ratio, as shown in Figure 4c. From Figure 3, it can be seen that in the ascending process of the intruder, the intruder actually experiences an alternating rise and fall. This is understandable as in each vibration cycle, the bed is lifted in the first half cycle, then the bottom wall moves downwards and so the particles also fall. Therefore, the filling of small particles in the void region underneath the big intruder is acknowledged as one of the major mechanisms for the overall ascending of the intruder [21], and the change of va with Ad is probably also related to this mechanism. The increase of va with increasing Ad has also been observed for a spherical intruder by Ahmad et al. [26] and by Rosato et al. [21]. As explained by Rosato et al., such an increase can be attributed to the fact that larger vibration amplitude From Figure 3, it can be seen that in the ascending process of the intruder, the intruder actually experiences an alternating rise and fall. This is understandable as in each vibration cycle, the bed is lifted in the first half cycle, then the bottom wall moves downwards and so the particles also fall. Therefore, the filling of small particles in the void region underneath the big intruder is acknowledged as one of the major mechanisms for the overall ascending of the intruder [21], and the change of v a with A d is probably also related to this mechanism. The increase of v a with increasing A d has also been observed for a spherical intruder by Ahmad et al. [26] and by Rosato et al. [21]. As explained by Rosato et al., such an increase can be attributed to the fact that larger vibration amplitude will cause more voids which promote the filling of small particles beneath the large intruder, and thus accelerate the ascending. This explanation is also applicable to the ascending of a non-spherical intruder.
However, when A > 1.5·d s , a further increase of vibration amplitude led to a decrease of v a . Such critical change can be linked with the structural changes in the vibrated bed. For spherical particles, Hsiao et al. [30] found that the segregation rate of a group of large intruders in a vibrated bed of small particles undergoes a nonmonotonic change with increasing vibration intensity. The primary cause for the critical change is that a too high vibration intensity can create very large voids, which let larger intruders fall down rather Processes 2021, 9, 448 7 of 15 than ascend. When vibration intensity is too high, the bed stops further expanding and the voids are so large that the probabilities for the larger intruder and smaller particles to fall down are equal, thus segregation no longer occurs. In addition, a later study [29] on the segregation of a single spherical intruder further proved that the ascending velocity has a strong dependence on the porosity variation of the particle bed.
Therefore, the non-monotonic response caused by A d for a non-spherical intruder can be linked to the variations of bed porosity (ε), which was investigated here. The porosity was calculated by dividing the particle bed into several horizontal layers with vertical thickness of 3.5·d s , and the porosity of each layer was calculated. The porosity variation (∆ε) is defined as the difference between the maximum and minimum of the bed porosity within one oscillation period, and the value is averaged over ten full-oscillation cycles. Figure 5 shows the porosity variation for different layers. Generally, ∆ε increased first and then decreased with the continuous increase of A d , and the critical changes occurred when A d = 1.5. The comparable change of v a to that of ∆ε hints that increasing the vibration amplitude would cause a larger ∆ε, which gives small particles more opportunities to fill voids beneath the intruder, thus accelerating the ascending. However, when vibration is too extensive, ∆ε reduces and the void filling becomes less efficient, and hence the ascending slows down.
cending of a non-spherical intruder.
However, when A > 1.5·ds, a further increase of vibration amplitude led to a decrease of va. Such critical change can be linked with the structural changes in the vibrated bed. For spherical particles, Hsiao et al. [30] found that the segregation rate of a group of large intruders in a vibrated bed of small particles undergoes a nonmonotonic change with increasing vibration intensity. The primary cause for the critical change is that a too high vibration intensity can create very large voids, which let larger intruders fall down rather than ascend. When vibration intensity is too high, the bed stops further expanding and the voids are so large that the probabilities for the larger intruder and smaller particles to fall down are equal, thus segregation no longer occurs. In addition, a later study [29] on the segregation of a single spherical intruder further proved that the ascending velocity has a strong dependence on the porosity variation of the particle bed.
Therefore, the non-monotonic response caused by Ad for a non-spherical intruder can be linked to the variations of bed porosity (ε), which was investigated here. The porosity was calculated by dividing the particle bed into several horizontal layers with vertical thickness of 3.5·ds, and the porosity of each layer was calculated. The porosity variation (∆ε) is defined as the difference between the maximum and minimum of the bed porosity within one oscillation period, and the value is averaged over ten full-oscillation cycles. Figure 5 shows the porosity variation for different layers. Generally, ∆ε increased first and then decreased with the continuous increase of Ad, and the critical changes occurred when Ad = 1.5. The comparable change of va to that of ∆ε hints that increasing the vibration amplitude would cause a larger ∆ε, which gives small particles more opportunities to fill voids beneath the intruder, thus accelerating the ascending. However, when vibration is too extensive, ∆ε reduces and the void filling becomes less efficient, and hence the ascending slows down. Similarly, the increase of va with the increase of the aspect ratio of the intruder can also be linked to the change of porosity variation at different bed heights. Figure 6a shows that ∆ε varied significantly with bed height, and higher layers showed smaller ∆ε. Therefore, when the centers of intruders with different aspect ratios are at the same height, as shown in Figure 6b, the elongated shape of the non-spherical intruder enables its lowerend to reach a deeper bed, where the void filling is more efficient. The more efficient void filling brings the non-spherical intruder of a larger aspect ratio a faster rise. Similarly, the increase of v a with the increase of the aspect ratio of the intruder can also be linked to the change of porosity variation at different bed heights. Figure 6a shows that ∆ε varied significantly with bed height, and higher layers showed smaller ∆ε. Therefore, when the centers of intruders with different aspect ratios are at the same height, as shown in Figure 6b, the elongated shape of the non-spherical intruder enables its lower-end to reach a deeper bed, where the void filling is more efficient. The more efficient void filling brings the non-spherical intruder of a larger aspect ratio a faster rise.

Effects of Vibration Frequency on v a
Similar to the effect of vibration amplitude, Figure 7 shows that, for both the spherical and the non-spherical intruders, the ascending velocity also increased and then decreased with increasing dimensionless vibration frequency. Comparably, non-spherical intruders have much higher v a than that of the spherical intruder. Here the dimensionless frequency is defined as f d = f ·(2·d s /g) 0.5 . Processes 2021, 9,

Effects of Vibration Frequency on va
Similar to the effect of vibration amplitude, Figure 7 shows that, for both the spherical and the non-spherical intruders, the ascending velocity also increased and then decreased with increasing dimensionless vibration frequency. Comparably, non-spherical intruders have much higher va than that of the spherical intruder. Here the dimensionless frequency is defined as fd = f‧(2·ds/g) 0.5 . Such a non-monotonic dependence was also observed in our previous research on binary mixtures, when A is fixed, and f varies [40]. The previous explanations of the influence of f were based on its impact on granular convection intensity [12,40]. However, as shown in Figure 8, there is no granular convection in the present research.

Effects of Vibration Frequency on va
Similar to the effect of vibration amplitude, Figure 7 shows that, for both the spherical and the non-spherical intruders, the ascending velocity also increased and then decreased with increasing dimensionless vibration frequency. Comparably, non-spherical intruders have much higher va than that of the spherical intruder. Here the dimensionless frequency is defined as fd = f‧(2·ds/g) 0.5 . Such a non-monotonic dependence was also observed in our previous research on binary mixtures, when A is fixed, and f varies [40]. The previous explanations of the influence of f were based on its impact on granular convection intensity [12,40]. However, as shown in Figure 8, there is no granular convection in the present research. Such a non-monotonic dependence was also observed in our previous research on binary mixtures, when A is fixed, and f varies [40]. The previous explanations of the influence of f were based on its impact on granular convection intensity [12,40]. However, as shown in Figure 8, there is no granular convection in the present research.
Similar to the above discussion on the effect of vibration amplitude, the non-monotonic change of v a may also be related to the change of bed porosity. As shown in Figure 9a, ∆ε of bed layer at a height of 20·d s increased first and then decreased with f d , which is consistent with the trend of v a . The change in ∆ε has a significant influence on filling interstices for small particles, and an increase or decrease of ∆ε would accordingly promote or suppress the ascending velocity. Besides, it is worth noting that the small particles can change from a disordered structure to an ordered structure under certain vibration conditions, which has been reported in our recent studies [41,42]. This can also impose an important influence on the intruder's ascending behavior and deserves further studies. ses 2021, 9, 448 9 of 16 Similar to the above discussion on the effect of vibration amplitude, the non-monotonic change of va may also be related to the change of bed porosity. As shown in Figure  9a, ∆ε of bed layer at a height of 20·ds increased first and then decreased with fd, which is consistent with the trend of va. The change in ∆ε has a significant influence on filling interstices for small particles, and an increase or decrease of ∆ε would accordingly promote or suppress the ascending velocity. Besides, it is worth noting that the small particles can change from a disordered structure to an ordered structure under certain vibration conditions, which has been reported in our recent studies [41,42]. This can also impose an important influence on the intruder's ascending behavior and deserves further studies.
In addition, according to [43], the influence of fd could also be linked to the change in the external energy input efficiency. The "granular temperature" that has been widely used to describe the kinetics of particle flow is linked with the random motion of particles [44,45]. Here the granular temperature (Tg) in the squared vicinity of the intruder was measured, and Tg is given by [46] T g = where Tgx, Tgy, and Tgz represent the granular temperature in x, y, and z directions, respectively; i refers to the ith particle and Ni is the total number of particles of interest; vix, viy, and viz denote the triaxial velocities of the ith particle and , , and are ensembleaverage velocities, respectively. Figure 9b shows that the variation of Tg with the change of fd is similar to that of va, suggesting a correlation between Tg and va. This indicates that an increase of fd will increase the random motion of small particles, which enhances the filling process and thereupon promotes the percolation. To the best of our knowledge, for the ascending of a single intruder in a vibrated bed, the non-monotonic dependence of the ascending velocity to vibration frequency has never been reported before. This may be due to the limited range of vibration frequency used in previous studies. For example, in the research by Khan Ahmad [26], only the cases of f > 50 Hz were investigated. Additionally, frequencies ranging in 4.07-4.99 Hz were investigated by Liffman et al. [22], who only found a positive effect of f on va when A is fixed. On the other hand, some studies only presented the effect of vibration frequency on va with fixed Γ [26,31], which simply showed a negative correlation of f on va.
These comparisons show that the effects of A and f on va are more complicated than previously reported, therefore they can be considered in a wider range and their individual influence can be more carefully examined. Thus, a phase diagram of va as a function of Ad and fd for our simulated results is presented as Figure 10, where the intruder is of dSau = 6 mm and a/b= 3.0. Note that the cases with extreme large Ad and fd were not included because the bed oscillates too violently to determine the position of the intruder. It can be seen that the phase diagram shows the symmetrical distribution of va with Ad and fd. Generally, a larger va can be obtained when Ad ∈ [1.25, 2.0] and fd < 0.6. The phase diagram also shows qualitative agreement with the variation trend of va against fd in [31], which reported that ascending velocity decreases with vibration frequency when Γ = 3.83 and Γ = 3.75. In addition, according to [43], the influence of f d could also be linked to the change in the external energy input efficiency. The "granular temperature" that has been widely used to describe the kinetics of particle flow is linked with the random motion of particles [44,45].
Here the granular temperature (T g ) in the squared vicinity of the intruder was measured, and T g is given by [46] T g = T gx + T gy + T gz 3 (7) where T gx , T gy , and T gz represent the granular temperature in x, y, and z directions, respectively; i refers to the ith particle and N i is the total number of particles of interest; v ix , v iy , and v iz denote the triaxial velocities of the ith particle and v x , v x , and v x are ensembleaverage velocities, respectively. Figure 9b shows that the variation of T g with the change of f d is similar to that of v a , suggesting a correlation between T g and v a . This indicates that an increase of f d will increase the random motion of small particles, which enhances the filling process and thereupon promotes the percolation.
To the best of our knowledge, for the ascending of a single intruder in a vibrated bed, the non-monotonic dependence of the ascending velocity to vibration frequency has never been reported before. This may be due to the limited range of vibration frequency used in previous studies. For example, in the research by Khan Ahmad [26], only the cases of f > 50 Hz were investigated. Additionally, frequencies ranging in 4.07-4.99 Hz were investigated by Liffman et al. [22], who only found a positive effect of f on v a when A is fixed. On the other hand, some studies only presented the effect of vibration frequency on v a with fixed Γ [26,31], which simply showed a negative correlation of f on v a .
These comparisons show that the effects of A and f on v a are more complicated than previously reported, therefore they can be considered in a wider range and their individual influence can be more carefully examined. Thus, a phase diagram of v a as a function of A d and f d for our simulated results is presented as Figure 10, where the intruder is of d Sau = 6 mm and a/b= 3.0. Note that the cases with extreme large A d and f d were not included because the bed oscillates too violently to determine the position of the intruder. It can be seen that the phase diagram shows the symmetrical distribution of v a with A d and f d . Generally, a larger v a can be obtained when A d ∈ [1.25, 2.0] and f d < 0.6. The phase diagram also shows qualitative agreement with the variation trend of v a against f d in [31], which reported that ascending velocity decreases with vibration frequency when Γ = 3.83 and Γ = 3.75.

Modeling of Ascending Velocity
It is shown in the above discussion that the inflection points in non-monotonic response curves (va-A and va-f) may vary with A and f simultaneously. Thus, before modeling, the critical change condition should be identified first. Several series of data with either A or f varied are selected to be analyzed together, which are listed in Table 4. Previously, Fernando et al. [29] found that for a spherical intruder, the critical change of va can be linked with the vibration velocity amplitude Vb, given by Vb = 2πAf(ds·g) −0.5 . In particular, va increases with Vb when Vb < 2.8, but decreases with Vb when Vb > 2.8. Similarly, the maximum va occurs when Vb = 2.75 for some of our series of results, as shown in Figure  11a. However, this is not always the case for other simulation results, as shown in Figure  11a. Instead, it is found that Γ = 4.83 can more uniformly identify the critical change of va for our simulation results, as shown in Figure 11b. Table 4. Parameters used in cases in Figure 11.

No.
fd A d dl/ds a/b Figure 10. v a as a function of f d and A d for the intruder of d Sau = 6 mm, a/b = 3.0. and • show the vibration conditions used in [31].

Modeling of Ascending Velocity
It is shown in the above discussion that the inflection points in non-monotonic response curves (v a -A and v a -f ) may vary with A and f simultaneously. Thus, before modeling, the critical change condition should be identified first. Several series of data with either A or f varied are selected to be analyzed together, which are listed in Table 4. Previously, Fernando et al. [29] found that for a spherical intruder, the critical change of v a can be linked with the vibration velocity amplitude V b , given by V b = 2πAf (d s ·g) −0.5 . In particular, v a increases with V b when V b < 2.8, but decreases with V b when V b > 2.8. Similarly, the maximum v a occurs when V b = 2.75 for some of our series of results, as shown in Figure 11a. However, this is not always the case for other simulation results, as shown in Figure 11a.
Instead, it is found that Γ = 4.83 can more uniformly identify the critical change of v a for our simulation results, as shown in Figure 11b. Table 4. Parameters used in cases in Figure 11.   Table 4.

No
Even the critical point at Γ = 4.83 is identified, and one can see that the data points do not collapse into one curve either in Figure 11a or b, indicating both Vb and Γ cannot be used as a single vibration parameter to model va. Comparing Vb and Γ, it can be noticed that they are both in proportion to Af n , but with different n values. Therefore, in our previous study, we have searched n for the best correlation between va and A d ⋅f d n in a binary system of spherical particles, and it was shown that A d ⋅f d 0.5 can be used to model va uniformly. Here, it is applied to our simulation results and other experiment results. Figure  12 shows that A d ⋅f d 0.5 can always be used to uniformly describe va for a set of data with different vibration conditions, though the dependency is different in different sets, which should be due to the differences in other parameters, such as different particles used [22,[31][32][33].  Table 4.
Even the critical point at Γ = 4.83 is identified, and one can see that the data points do not collapse into one curve either in Figure 11a,b, indicating both V b and Γ cannot be used as a single vibration parameter to model v a . Comparing V b and Γ, it can be noticed that they are both in proportion to Af n , but with different n values. Therefore, in our previous study, we have searched n for the best correlation between v a and A d · f n d in a binary system of spherical particles, and it was shown that A d · f 0.5 d can be used to model v a uniformly. Here, it is applied to our simulation results and other experiment results. Figure 12 shows that A d · f 0.5 d can always be used to uniformly describe v a for a set of data with different vibration conditions, though the dependency is different in different sets, which should be due to the differences in other parameters, such as different particles used [22,[31][32][33].
In addition, Figure 13a shows that the shape of the intruder has a significant impact on its ascending. This seems to be different from the conclusion in [26] that the shape of the intruder has little effect. The reason is probably because that in the experiments in [26] there was strong convection in the vibrated bed, but in our simulations with PBCs the convection was not formed. In such conditions, the shape of the intruder should be considered. In fact, as shown in Figure 13a, v a is proportional to Φ γ , where Φ is the sphericity of the intruder and γ is the exponent. formly. Here, it is applied to our simulation results and other experiment results. Figure  12 shows that A d ⋅f d 0.5 can always be used to uniformly describe va for a set of data with different vibration conditions, though the dependency is different in different sets, which should be due to the differences in other parameters, such as different particles used [22,[31][32][33]. In addition, Figure 13a shows that the shape of the intruder has a significant impact on its ascending. This seems to be different from the conclusion in [26] that the shape of the intruder has little effect. The reason is probably because that in the experiments in [26] there was strong convection in the vibrated bed, but in our simulations with PBCs the convection was not formed. In such conditions, the shape of the intruder should be considered. In fact, as shown in Figure 13a, va is proportional to Φ γ , where Φ is the sphericity of the intruder and γ is the exponent. Based on the above discussion, a model is proposed to predict the ascending velocity of a non-spherical intruder in a vibrated bed of small mono-size spherical particles: where K1 and K2 are proportionality coefficients, ε1 and ε2 are correction coefficients, respectively. These coefficients can be fitted for a given granular system. As shown in Figure  13b, the better fitting of our simulation results can be obtained when K1, K2, ε1, and ε2 are 0.095, −0.04, −0.1, and 0.23, respectively. As the ascending velocity may be changed when using different equivalent diameters, here γ for using dSau, dArea, and dVol are respectively given as γSau = −6.74, γArea = −6.01, and γVol = −6.31. It can be seen that there is a slight difference between γSau, γVol, and γArea, and va calculated with the same sphericity but different equivalent diameters have the highest value when dSau is used, while the lowest value when dArea is used. To further test the model, the experimental results for a single spherical intruder (Φ = 1.0) in [31][32][33] were also fitted by our Equation (9). The predicted values by the fitting equations were compared to the experimental data in Figure 14. The predictions are in good agreement with the experimental values, showing the model has general applicability. Based on the above discussion, a model is proposed to predict the ascending velocity of a non-spherical intruder in a vibrated bed of small mono-size spherical particles: where K 1 and K 2 are proportionality coefficients, ε 1 and ε 2 are correction coefficients, respectively. These coefficients can be fitted for a given granular system. As shown in Figure 13b, the better fitting of our simulation results can be obtained when K 1 , K 2 , ε 1 , and ε 2 are 0.095, −0.04, −0.1, and 0.23, respectively. As the ascending velocity may be changed when using different equivalent diameters, here γ for using d Sau , d Area , and d Vol are respectively given as γ Sau = −6.74, γ Area = −6.01, and γ Vol = −6.31. It can be seen that there is a slight difference between γ Sau , γ Vol , and γ Area , and v a calculated with the same sphericity but different equivalent diameters have the highest value when d Sau is used, while the lowest value when d Area is used.
To further test the model, the experimental results for a single spherical intruder (Φ = 1.0) in [31][32][33] were also fitted by our Equation (9). The predicted values by the fitting equations were compared to the experimental data in Figure 14. The predictions are in good agreement with the experimental values, showing the model has general applicability.

Conclusions
In this paper, the segregation of a large intruder of different shapes in a vibrated bed of small particles was investigated by a DEM model. Several important conclusions can be drawn from the results.
Firstly, the ascending velocity (va) of a non-spherical intruder increases first and then decreases with the rising of the dimensionless vibration amplitude (Ad). There is a similar variation trend of va with increasing dimensionless vibration frequency (fd). An intruder of larger aspect ratio rises faster under the same vibration conditions.
Secondly, under both the influences of Ad and fd, the change of the porosity variation of the vibrated bed during a vibration period (∆ε) is similar to that of va. This indicates that the changes of va can be primarily attributed to the change in void filling of small particles, as in our simulations the convection was not generated, and the percolation mechanism dominated the segregation. This can also be used to explain the effect of intruder shape, as a more elongated intruder can make void filling more efficient, resulting in a larger va. The change of granular temperature (Tg) also shows a link to the change of va.
Thirdly, not only for current simulations but also for the experimental data in the literature, A d ⋅f d 0.5 is shown to be an appropriate combined parameter to characterize the influence of A and f on va. Additionally, va is related to the negative nth power of the intruder's sphericity Φ. Based on these findings, a piecewise equation has been proposed to model va as a function of Ad, fd, and Φ. The equation can be well fitted with all our simulation results and the experimental results in the literature. The coefficients in the equation are dependent on other particle properties, which deserve further study.

Conclusions
In this paper, the segregation of a large intruder of different shapes in a vibrated bed of small particles was investigated by a DEM model. Several important conclusions can be drawn from the results.
Firstly, the ascending velocity (v a ) of a non-spherical intruder increases first and then decreases with the rising of the dimensionless vibration amplitude (A d ). There is a similar variation trend of v a with increasing dimensionless vibration frequency (f d ). An intruder of larger aspect ratio rises faster under the same vibration conditions.
Secondly, under both the influences of A d and f d , the change of the porosity variation of the vibrated bed during a vibration period (∆ε) is similar to that of v a . This indicates that the changes of v a can be primarily attributed to the change in void filling of small particles, as in our simulations the convection was not generated, and the percolation mechanism dominated the segregation. This can also be used to explain the effect of intruder shape, as a more elongated intruder can make void filling more efficient, resulting in a larger v a . The change of granular temperature (T g ) also shows a link to the change of v a .
Thirdly, not only for current simulations but also for the experimental data in the literature, A d · f 0.5 d is shown to be an appropriate combined parameter to characterize the influence of A and f on v a . Additionally, v a is related to the negative nth power of the intruder's sphericity Φ. Based on these findings, a piecewise equation has been proposed to model v a as a function of A d , f d , and Φ. The equation can be well fitted with all our simulation results and the experimental results in the literature. The coefficients in the equation are dependent on other particle properties, which deserve further study. Funding: This research was funded by the joint Ph.D. program of "double first rate" construction disciplines of CUMT.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to the protection of ongoing original research.