separations

Article
The Nonlinear Dynamic Behavior of a Particle on a Vibrating
Screen Based on the Elastoplastic Contact Model
Deyi He, Chusheng Liu * and Sai Li

School of Mechatronic Engineering, China University of Mining and Technology, Xuzhou 221008, China
* Correspondence: [email protected]

Abstract: The dynamic response of particles is closely related to screening efficiency. To study
the dynamic response of particles, the dynamic equations of a particle on a screening surface are
established based on the elastoplastic contact model of spherical particles and are solved for the coal
particle. Then, the trajectories of the particles are given with different falling heights and particle
radii. The completely different trajectories with slight changes in the falling height and particle radius
indicate strong nonlinearity. Second, the nonlinear dynamic behavior under different amplitudes and
frequencies is discussed, and the route of transition from quasiperiodic motion to chaotic motion is
revealed. Finally, we discuss the average speed along the screening surface considering the frequency,
amplitude, friction coefficient, inclination angle, and vibration direction angle. In addition, the
convergence conditions of particle motion are proposed, and they are only affected by the inclination
angle and friction angle. The results show that in the normal direction of the vibrating screen, the
particle motion is quasiperiodic at low frequencies. With increasing frequency, the motion of the
particle becomes chaotic, and its Poincaré map becomes petal-shaped. In addition, the number of
petals increases at the mutation of the bifurcation diagram. The increase in frequency, amplitude and
inclination angle and the decrease in friction coefficient lead to an increase in particle speed along the
screen surface. In addition, the particle speed reaches a maximum when the vibration direction angle
is 65◦ . This work provides a theoretical basis for controlling the thickness of granular material flow
Citation: He, D.; Liu, C.; Li, S. The on a vibrating screen and selecting screening process parameters.
Nonlinear Dynamic Behavior of a
Particle on a Vibrating Screen Based Keywords: Hertz contact; nonlinear dynamic; vibrating screen; chaos motion; granular coal
on the Elastoplastic Contact Model.
Separations 2022, 9, 216. https://
doi.org/10.3390/separations9080216

Academic Editor: Qicheng Feng 1. Introduction

High particle efficiency is the most important problem for vibrating screens [1]. The
Received: 23 July 2022
present studies for the screening theory of granular materials mainly focus on simulating the
Accepted: 9 August 2022
screening process by using the discrete element method (DEM) [2], the screening experiment
Published: 12 August 2022
for choosing optimal vibrating parameters (amplitude, frequency, inclination angle, etc.) [3],
Publisher’s Note: MDPI stays neutral the trajectory of particles bouncing on a vibrating screen and the invention of new types of
with regard to jurisdictional claims in vibrating screens [4]. The ultimate purpose is to improve screen efficiency.
published maps and institutional affil- The essence of screening is to apply a small amplitude and suitable vibrational fre-
iations.
quency to make the particles bounce and filter smaller particles under the screen [5].
Therefore, the key problem of the screening process is the dynamic behavior of the par-
ticles on the vibrating screen. DEM, as an effective numerical simulation technique, is a
powerful method for predicting the motion and statement of the particles in the screening
Copyright: © 2022 by the authors.
Licensee MDPI, Basel, Switzerland.
process [6,7]. DEM can obtain the key parameters, which include screening, particle speed,
This article is an open access article
mass distribution, and impact force. The screening mechanism is revealed more easily
distributed under the terms and with the help of DEM simulations. EDEM (the world’s first multi-purpose DEM modeling
conditions of the Creative Commons software, which can be used to simulate and analyze the production process of particle pro-
Attribution (CC BY) license (https:// cessing and manufacturing equipment in industrial production), the most commonly used
creativecommons.org/licenses/by/ simulation software for DEM, is a reliable method for simulating the screening processing
4.0/). of vibrating screening.

Separations 2022, 9, 216. https://doi.org/10.3390/separations9080216 https://www.mdpi.com/journal/separations

Separations 2022, 9, 216 2 of 14

However, the dynamic behavior of the particles in EDEM software does not strictly obey
the law of Hertz’s contact model and Thornton’s elastoplastic contact model. The coefficient
of restitution in EDEM is set as a constant, which ignores the nonlinearity of the contact model
despite the calculation of the contact force being based on Hertz’s model because of the large
number of particles. Furthermore, the calculation of the contact force is not important for
the screening process. In addition, the nonlinear dynamic behavior of particles on the screen
surface reflects the screening process under different screening parameters and granular
materials, and it is of great significance for revealing the screening mechanism.
The experimental method is a powerful way to study the screening mechanism as well,
and it mainly focuses on screening granular materials to measure the content on and under
the screen to calculate the screening efficiency. Meanwhile, the screening parameters are
optimized by repeatedly changing the vibration parameters. Jiang [8] used a disequilibrium
excitation vibrating screen to screen moist coal to improve the screening efficiency by
changing the vibrating parameters. The physical screening method is absolutely accurate,
but it is labor-intensive and time-consuming. The adjustable parameters are limited by the
equipment. Additional equipment, including a high-speed camera and three-component
acceleration sensors, should be set when extra parameters, such as particle speed and
vibrating trajectory, are needed. In addition, the small particles are difficult to monitor
by using a high-speed camera because the particle flow covers the small particles in the
screening process.
Therefore, the behavior of particles on the screening surface is addressed by using a
theoretical method. Wang [9] discussed one particle’s dynamic behavior on a nonsimple
harmonic vibrating screen, and the particle speed in different areas of the vibrating screen
was given. However, the contact model between the vibrating screen and particles was linear.
In the field of vibration screening, little research considers the nonlinear contact model
between vibrating screens and granular materials. Yang [10] discussed the vibroimpact
response of the particle elastic impact on a metal plate without vibration based on Hertz’s
contact model. Dong [11] discussed the dynamic behavior of a vibrating plate at high
frequency by extracting the particle load based on a linear contact model.
The nonlinear contact model is the fundamental reason for the nonlinear motion of
particles on a vibrating screen. In 1896, Hertz [12] derived the elastic contact model between
two spheres. The force–displacement relationship (FDR) is given by

4 √ 1.5
F= E Rδ (1)
3
where E is Young’s modulus, R is the radii of particles, δ is relative deformation.
Hurt and Crossley [13] proposed the viscoelastic contact model of two spheres based
on Hertz’s theory, and the FDR is shown as follows:
.
F = Kδn + χδn δ (2)

where K is stiffness coefficient, χ is damping coefficient, n is the coefficient of contact

according to the contact type.
In 1997, Thornton [14] found that energy loss caused by elastoplastic deformation
is the main reason for different coefficients of restitution between different spheres. The
collision process can be divided into three periods, including the elastic compression period,
elastoplastic compression period, and elastic recovery period, and the FDR is given by:

4 √
F = πδy R(δ − δy ) + E Rδy1.5 (3)
3
4 √
1.5
F= E R δ − δy (4)
3
where δy is the ultimate deformation in the elastic contact period.
 E R ( −  y )
4 1.5
F= (4)
3
Separations 2022, 9, 216 where δy is the ultimate deformation in the elastic contact period. 3 of 14
Many researchers have focused on the theoretical model of particle contact based on
the above theory [15–23]. Essentially, there are few differences between these calculations
with different models [24]. Meanwhile, the validity of the theoretical models is difficult to
Many researchers have focused on the theoretical model of particle contact based on
prove because of the lack of experiments.
the above theory [15–23]. Essentially, there are few differences between these calculations
In this research, the dynamic equations of a particle bouncing on a linear vibrating
with different models [24]. Meanwhile, the validity of the theoretical models is difficult to
screen surface are established based on Thornton’s contact model. We chose coal particles,
prove because of the lack of experiments.
which are the most widely granular material being screened, as screening material. The
In this research, the dynamic equations of a particle bouncing on a linear vibrating
screening
screen surface
surface is treated as
are established a rigid
based body because
on Thornton’s of the
contact large We
model. differences between
chose coal the
particles,
iron screen surface and coal particles in terms of the elastic modulus. Eventually,
which are the most widely granular material being screened, as screening material. The we give
the nonlinear
screening behavior
surface of a as
is treated particle
a rigidinbody
different directions.
because of the large differences between the
iron screen surface and coal particles in terms of the elastic modulus. Eventually, we give
2. Dynamic
the nonlinear Equations
behavior ofofa aparticle
Particleinon Vibrating
different Screening
directions.
The dynamic model of the particle-vibrating screening (P–VS) system is shown in
2.Figure
Dynamic Equations
1. The of a Particle
particle bounces on theonsinusoidal
Vibratingvibration
Screening screening surface. Furthermore,
the The
movement
dynamicconsists
modelofoffour
theperiods:
particle-vibrating screening (P–VS) system is shown in
Figure
1. The1. The particle
particle hasbounces
no contacton with
the sinusoidal
the screenvibration
surface andscreening surface. Furthermore,
flies freely.
the
2. movement consists
The particle of four
impacts periods:
the screen surface with elastic loading.
1.3. The
Theparticle
particlehas
impacts the screen
no contact with surface with
the screen elastoplastic
surface and fliesloading.
freely.
2.4. The
Theparticle
particleimpacts
impactsthe
thescreen
screensurface
surfacewith
withelastic
elasticloading.
unloading and returns to period
3. 1. particle impacts the screen surface with elastoplastic loading.
The
4. The particle impacts the screen surface with elastic unloading and returns to period 1.

F
a

Y fcF
Asin(ωt)
X

Screening g
surface
θ

Figure1.1.Dynamic
Figure Dynamicmodel
modelofofthe
theparticle-vibrating
particle-vibratingscreening
screeningsystem.
system.
On the basis of Thornton and Hertz’s theory, we establish the dynamic equations of
On the basis of Thornton and Hertz’s theory, we establish the dynamic equations of
particles on vibrating screens by considering the elastoplastic contact between granular
particles on vibrating screens by considering the elastoplastic contact between granular
materials and vibrating screens. For these purposes, the following assumptions are made:
materials and vibrating screens. For these purposes, the following assumptions are made:
1. The screen surface of the vibrating screen is rigid.
1.
2. Thetangential
The screen surface of the vibrating
displacement screenmaterials
of granular is rigid. is ignored, and only the fraction
2. isThe tangential displacement of granular materials is ignored, and only the fraction is
considered.
3. considered.
The rolling of granular materials is ignored.
3. The rolling of granular materials is ignored.
Therefore, the equivalent radius R and Young’s modulus E can be given as
Therefore, the equivalent radius R and Young’s modulus E can be given as
E = E1 , R = R1 (5)
E = E1 , R = R1 (5)
where R1 is the radius of the particle, and E1 is Young’s modulus of the particle. According
to the above theories, the dynamic equations can be written as
..
my = −mg · cos θ + F (yv − y) (6)
.. . .
m x = mg · sin θ + sgn( xV − x ) f c F (yv − y) (7)
where yv − y is the relative deformation of the particle and screen surface and only occurs at
contact time, yv and xv is the displacement of the vibrating screen in normal and tangential
Separations 2022, 9, 216 4 of 14

direction, y and x is the displacement of the particle in normal and tangential direction, m
is the mass of the particle, θ is the inclination angle of the vibrating screen, and fc is the
friction coefficient. Hence, the contact force F(yv − y) can be written as Equations (8)–(11)
based on Thornton’s work
F (yv − y) = 0, y ≥ yv (8)
4 √
F (yv − y) = E R(yv − y)1.5 , y < yv &yv − y ≤ δy (9)
3
. .
F (yv − y) = πδy R(yv − y − δy ) + Py , y < yv &yv − y ≤ δy &yv − y < 0 (10)
4 √
1.5 . .
F (yv − y) = E R yv − y − δy , y < y v &y v − y > 0 (11)
3
where Py is the elastic contact force.

π 3 σy3 R2
Py = (12)
6E

π 2 σy2 R
δy = (13)
4E2
where δy is the maximum deformation in the elastic load period, and σy is the ultimate
pressure, which is defined as 2.5 times the compressive strength of the granular material [25]

σy = 2.5σ (14)

In a normal vibrating screen, the kinematic equations can be defined as


xv = A sin(ωt) · sin a
(15)
yv = A sin(ωt) · cos a

where A is the amplitude of the vibrating screen, ω is the rotational speed, and a is the
inclination angle. Using Equations (6)–(11) and (15), the dynamic response can be obtained.

3. Results and Discussion

3.1. Material Properties of the P–VS System
The coal particle is chosen for solving the dynamic equations, and the particle mechan-
ical properties and vibrating parameters are shown in Table 1.

Table 1. Particle mechanical properties and vibrating parameters.

Property Value
Density (kg/m3 ) 1300
Elastic modulus (Pa) 3.5 × 109
Compressive stress (Mpa) 30
Radius (mm) 3
Frequency of vibrating screen f (Hz) 8–26
Amplitude of vibrating screen A (mm) 1–10
Inclination angle of vibrating screen θ (◦ ) 6–30
Vibration mode Linear
Friction coefficient fc 0.5

3.2. The Dynamic Behavior with Different Falling Heights and Particle Radii
The dynamic behavior in the normal direction is independent of that in the tangential
direction but not vice versa through Equations (6) and (7). Hence, the particle dynamic
behavior of a circular vibrating screen and linear vibrating screen in the normal direction is
the same with the same amplitude in the normal direction. Figure 2a shows the particles
dropping from three heights (0.3, 0.301, and 0.302 m) under f = 16 Hz, and A = 4 × 10−3 m.
 3.2. The Dynamic Behavior with Different Falling Heights and Particle Radii
The dynamic behavior in the normal direction is independent of that in the tangential
direction but not vice versa through Equations (6) and (7). Hence, the particle dynamic
behavior of a circular vibrating screen and linear vibrating screen in the normal direction
Separations 2022, 9, 216
is the same with the same amplitude in the normal direction. Figure 2a shows the particles 5 of 14
dropping from three heights (0.3, 0.301, and 0.302 m) under f = 16 Hz, and A = 4 × 10−3 m.
Figure 2b shows the frequency spectra of the dynamic responses of the particle. The trends
Figurethe
are basically 2b same,
showsbut
the the
frequency spectra
curves of of the dynamic
the frequency spectraresponses of the particle.
at low frequencies The trends
are quite
are As
different. basically the same,
the figure shows,but thethe curvesmotion
particle of the frequency spectra at low
becomes completely frequencies
different are quite
after 2–3
bounces,different. As the figure
which implicates the shows,
strong the particlephenomenon
nonlinear motion becomes of thecompletely different after
particle-vibrating
screen 2–3 bounces,
system. Smallwhich implicates
changes in heightthe strong
lead nonlinear
to great changesphenomenon of the particle-vibrating
in particle trajectories, which
screen system. Small changes
is consistent with actual particle motion. in height lead to great changes in particle trajectories, which
is consistent with actual particle motion.

(a)

(b)
Figure 2. Dynamic
Figure response
2. Dynamic of the particles
response with different
of the particles falling heights:
with different (a) the trajectories
falling heights: of the of the
(a) the trajectories
particles and (b)and
particles the frequency spectraspectra
(b) the frequency of the particles.
of the particles.

In addition, the radius of the particles affects the trajectory. As shown in Figure 3, the
trajectories for three particle radii are given with the above parameters. The trajectories
of the three particles begin to differ after seven bounces. These results prove that minor
changes in any parameter might lead to large differences in the P–VS system. In the past,
the dynamic equations of vibrating screens only considered the excitation force given by
motors. Based on our work, the impact force of particles can be considered in the dynamic
equations because the accurate impact time is given for particles with different ranges of
radii and falling heights.
 the three particles begin to differ after seven bounces. These results prove that minor
changes in any parameter might lead to large differences in the P–VS system. In the past,
the dynamic equations of vibrating screens only considered the excitation force given by
motors. Based on our work, the impact force of particles can be considered in the dynamic
Separations 2022, 9, 216 equations because the accurate impact time is given for particles with different ranges 6ofof 14
radii and falling heights.

Figure 3. Dynamic
Figure responses
3. Dynamic of particles
responses with
of particles different
with radii.
different radii.

3.3.3.3. Formatting
Formatting of Mathematical
of Mathematical Components
Components
TheThe amplitude
amplitude of most
of most vibrating
vibrating screens
screens is between
is between 1–101–10
mm. mm. Therefore,
Therefore, thethe dynamic
dynamic
behavior at 1–10 mm amplitude and 16 Hz frequency is discussed,
behavior at 1–10 mm amplitude and 16 Hz frequency is discussed, in which 16 Hz is the in which 16 Hz is the
common frequency of normal
common frequency of normal vibrating vibrating screens. Figure 4a shows the bifurcation diagram of
screens. Figure 4a shows the bifurcation diagram
different amplitudes at θ = 10 ◦ , a = 90◦ , f = 16 Hz, and f = 0.5. With increasing amplitude,
of different amplitudes at θ = 10°, a = 90°, f = 16 Hz, and fc = c0.5. With increasing amplitude,
thethe dynamic
dynamic system
system of P–VS
of P–VS becomes
becomes moremore unstable.
unstable. Moreover,
Moreover, thethe attraction
attraction domain
domain
becomes larger. In Figure 4a–i, the phase diagrams and the Poincaré maps under different
becomes larger. In Figure 4a–i, the phase diagrams and the Poincaré maps under different
amplitudes are shown. When the amplitude is 1 mm, the phase diagram shows that the
amplitudes are shown. When the amplitude is 1 mm, the phase diagram shows that the
particles bounce once on the screening surface in each period, and the collision occurs when
particles bounce once on the screening surface in each period, and the collision occurs
the particle position is 0.0003–0.0009 m. The points on the Poincaré map are drawn around
when the particle position is 0.0003–0.0009 m. The points on the Poincaré map are drawn
the particle velocity of −0.17 m/s and form a pattern similar to a nebula, indicating the
around the particle velocity of −0.17 m/s and form a pattern similar to a nebula, indicating
chaotic motion of the system. When the amplitude is 1.5 mm, the chaotic attractor changes
the chaotic motion of the system. When the amplitude is 1.5 mm, the chaotic attractor
greatly. The bouncing times between the particle and the screen surface decrease in one
changes greatly. The bouncing times between the particle and the screen surface decrease
period, which may occur at any time within the period of vibrating screen movement. In
in one period, which may occur at any time within the period of vibrating screen move-
the phase diagram, the hole in the middle of the chaotic attractor disappears. In addition,
ment. In the phase diagram, the hole in the middle of the chaotic attractor disappears. In
the Poincaré map changes from nebular to half-flower shaped with three petals, and the
addition, the Poincaré map changes from nebular to half-flower shaped with three petals,
value range of the Poincaré map is wider. Increasing the amplitude from 1.5 mm to 5 mm
andhasthealmost
value range of the
no effect on Poincaré
the attractor mapofis the
wider. Increasing
system, and the thePoincaré
amplitude mapfrom 1.5 mm
changes from
to 5three
mm has almost no effect on the attractor of the system, and the
petals to seven petals. When the amplitude increases to 9 mm, the Poincaré map Poincaré map changes
from three to
changes petals to seven
12 petals. petals. When
In conclusion, when thetheamplitude
amplitude increases
of the P–VSto 9 system
mm, the is Poincaré
1.5 mm, the
map changes to 12 petals. In conclusion, when the amplitude of the
chaotic attractor of the system is small, and the shape of the Poincaré map is nebulous. P–VS system is 1.5
With
mm, the chaotic attractor of the system is small, and the shape of the
increasing amplitude, the motion range of the particles increases, and the chaotic attractor Poincaré map is
nebulous. With increasing amplitude, the motion range of the particles
of the system increases. At the same time, the Poincaré map becomes petal-shaped, and the increases, and the
chaotic
largerattractor of the is,
the amplitude system
the moreincreases. At the
petals there aresame
on thetime, the Poincaré
Poincaré map. When map thebecomes
amplitude
petal-shaped, and the larger the amplitude is, the more petals
reaches 10 mm, the number of petals reaches 13 and tables should be cited in the there are on the Poincaré
main text
map. When 1,
as Figure theTable
amplitude
1, etc. reaches 10 mm, the number of petals reaches 13 and tables
should be cited in the main text as Figure 1, Table 1, etc.
 Separations 2022, 9, x FOR PEER REVIEW 7 of 15
Separations 2022, 9, 216 7 of 14

(a)

(b) (c)

(d) (e)

Figure 4. Cont.
 Separations 2022, 9, x FOR PEER REVIEW 8 of 15
Separations 2022, 9, 216 8 of 14

(f) (g)

(h) (i)
Figure4.4.Nonlinear
Figure Nonlinear behavior
behavior of of a particle
a particle at different
at different amplitudes:
amplitudes: (a) (a) Bifurcation
Bifurcation diagram
diagram of differ-
of different
ent amplitudes; (b) Phase diagram at A = 1 mm; (c) Poincaré map at A = 1 mm; (d) Phase diagram at
amplitudes; (b) Phase diagram at A = 1 mm; (c) Poincaré map at A = 1 mm; (d) Phase diagram at
A = 1.5 mm; (e) Poincaré map at A = 1.5 mm; (f) Phase diagram at A = 5 mm; (g) Poincaré map at A
A = 1.5 mm; (e) Poincaré map at A = 1.5 mm; (f) Phase diagram at A = 5 mm; (g) Poincaré map at
= 5 mm; (h) Phase diagram at A = 9 mm; (i) Poincaré map at A = 9 mm.
A = 5 mm; (h) Phase diagram at A = 9 mm; (i) Poincaré map at A = 9 mm.
3.4.Effects
3.4. EffectsofofFrequency
Frequencyononthe theDynamic
DynamicBehavior
BehaviorofofP–VS
P–VSSystem
System
Thedynamic
The dynamicbehavior behavior at at 44 mm
mm amplitude
amplitudeand and6–266–26Hz Hzfrequency
frequency is is
discussed
discussedin this
in
section. Figure 5a shows the bifurcation diagram of different
this section. Figure 5a shows the bifurcation diagram of different frequencies at A = 4 mm, frequencies at A = 4 mm, θ=
10°, ◦
a = 90°, ◦
and f = 0.5. According to Figure 5b,c, the particle motion
θ = 10 , a = 90 , and fc = 0.5. According to Figure 5b,c, the particle motion exhibits quasiperi-
c exhibits quasiperiodic
motion
odic motionat aatfrequency
a frequency of 6ofHz,6 Hz,andand
thetheparticle bounces
particle bounces 5 times each
5 times period,
each period,in which
in whichone
oneperiod
periodis is
approximately
approximately1/6 1/6s. s.Moreover,
Moreover,the thePoincaré
Poincaré mapmap shows aa circular
circulardot dotarea
area
withiny =
within y =0.95
0.95−−1.17 1.17××10 10−−33 m.
m. When
When the frequency increasesincreases to to6.5
6.5Hz,
Hz,thetheP–VS
P–VSsystem
system
changesfrom
changes fromquasiperiodic
quasiperiodicmotion motiontotochaotic
chaoticmotion
motionininwhichwhichbifurcation
bifurcationoccurs.
occurs.As As
shown in Figure 5d,e, the particle bounces on the screening
shown in Figure 5d,e, the particle bounces on the screening surface at any moment, and surface at any moment, and
itsitsPoincaré
Poincarémap mapresembles
resemblesa asemicircular
semicircularring ringabove
abovea asmall
smallfoot.
foot.When
Whenthe thefrequency
frequency
increasestoto88Hz,
increases Hz, its its nonlinear
nonlinear dynamic behavior behavior isisthe thesame
sameasasthe theabove
above content
contentof of
the
theamplitude
amplitude (Section
(Section 3.2), in in
3.2), which
which thethePoincaré
Poincaré map
map shows
shows a apetal
petalflower,
flower,asasshown
shownin
inFigure
Figure5f,g.
5f,g.InInaddition,
addition, each
each mutation
mutation ononthethe bifurcation
bifurcation diagram
diagram leadsleads to one
to one more more
petal
petal
on theon the Poincaré
Poincaré map.map. In In short,the
short, theP–VS
P–VSsystem
system presents quasiperiodic
quasiperiodicmotion motionwhenwhen
working
workingatatlow lowfrequencies.
frequencies.When Whenthe thefrequency
frequencyincreases
increasestoto6.5 6.5Hz,
Hz,the thebifurcation
bifurcationphe-phe-
nomenon
nomenonoccurs. occurs.At Atthis
thistime,
time,the thesystem
systempresents
presentschaotic
chaoticmotion,
motion,the theattractor
attractorbecomes
becomes
larger,
larger,andand thethelawlaw is is
difficult
difficult to to
find. With
find. Witha further
a furtherincrease in frequency,
increase in frequency, the attractor of
the attractor
theof system
the systembecomesbecomes semicircular,
semicircular, and itsandPoincaré map becomes
its Poincaré map becomes petal-shaped. After After
petal-shaped. this
stage, an increase in frequency will only increase the radius of the semicircle of the chaotic
attractor and the number of petals on the Poincaré map.
 Separations 2022, 9, x FOR PEER REVIEW 9 of 15

Separations 2022, 9, 216 this stage, an increase in frequency will only increase the radius of the semicircle9of the
of 14
chaotic attractor and the number of petals on the Poincaré map.

(a)

(b) (c)

(d) (e)

Figure 5. Cont.
Separations 2022, 9, x FOR PEER REVIEW 10 of 15
Separations 2022, 9, x FOR PEER REVIEW 10 of 15
Separations 2022, 9, 216 10 of 14

(f) (g)
(f) (g)Bifurcation diagram of differ-
Figure 5. Nonlinear behavior of a particle at different frequencies: (a)
ent frequencies;
Figure (b) Phase diagram at f = 6 Hz; (c) Poincaré map at f (a)= 6 Bifurcation
Hz; (d) Phase diagram at f =
Figure 5.5.Nonlinear
Nonlinear behavior
behavior ofofa aparticle
particleatat different
different frequencies:
frequencies: (a) Bifurcation diagram
diagram of differ-
of different
6.5 Hz; (e) Poincaré map at f = 6.5 Hz; (f) Phase diagram at f = 8 Hz; (g) Poincaré map at
ent frequencies; (b) Phase diagram at f = 6 Hz; (c) Poincaré map at f = 6 Hz; (d) Phase diagram at f = f = 8 Hz.
frequencies; (b) Phase diagram at f = 6 Hz; (c) Poincaré map at f = 6 Hz; (d) Phase diagram at
6.5 Hz; (e) Poincaré map at f = 6.5 Hz; (f) Phase diagram at f = 8 Hz; (g) Poincaré map at f = 8 Hz.
f = 6.5
3.5. TheHz; (e) Poincaré
Dynamic map of
Behavior at af =Particle
6.5 Hz;along
(f) Phase
the diagram
Screening f = 8 Hz; (g) Poincaré map at f = 8 Hz.
at Surface
3.5.The
3.5. Thethis
In Dynamic
section,
Dynamic Behavior ofofaaParticle
the particle
Behavior Particle along
speedalong
under the
the
the Screening
influence
Screening Surface
of amplitude, frequency, incli-
Surface
nation Inangle,
Inthis and
this section,friction
section, the coefficient
the particle
particle speed along
speed under the
underthescreening
theinfluence surface is mainly frequency,
influenceofofamplitude,
amplitude, discussed.
frequency, Fig-
incli-
in-
ure 6 shows
nation
clination angle, the
angle, anddisplacement
andfriction
friction of the along
coefficient
coefficientparticle
theat
along f =screening
screening
the 16 Hzsurface
and Ais= mainly
surface 4ismm. Basically,
discussed.
mainly the
Fig-
discussed.
displacement
ure 6 shows
Figure 6 shows of
the the
the particle keeps
displacement
displacement increasing
ofofthe
theparticle
particleatatataffconstant
==16 Hz speed,
16 Hz and
and A =in4which the average
mm. Basically,
mm. Basically, the
the
speed V
displacementis approximately
of the particle0.19 m/s.
keeps In addition,
increasing at each
a bounce
constant reduces
speed,
displacement of the particle keeps increasing at a constant speed, in which the average
a in the
which speed
the in the
average
tangential
speedVVaa is
speed direction, as shown
isapproximately
approximately inm/s.
0.19
0.19 Equations (6) and each
m/s. In addition,
addition, (7).
each bounce
bounce reduces
reduces the
the speed
speed in
inthethe
tangentialdirection,
tangential direction,asasshown
shownininEquations
Equations(6) (6)and
and(7). (7).

Figure
Figure 6.
6. Displacement
Displacementof
of the
the particle
particle at
at ff ==16
16Hz
Hzand
andAA==44mm.
mm.
Figure 6. Displacement of the particle at f = 16 Hz and A = 4 mm.
Figure 77 shows
Figure shows the average speed with with AA == 1–10 mm, f == 8–26 8–26 Hz, 10◦ , and
Hz, θθ == 10°, 90◦ .
and aa ==90°.
In general,
In general,
Figure the
the average
average
7 shows the speed
speed
averageof of
the the particle
particle
speed A increases
withincreases
= 1–10 with
mm, with increasing
fincreasing
= 8–26 θ =amplitude
Hz,amplitude
10°, andand and
a =fre-
90°.
frequency.
quency. At At
In general, low
thelow amplitude
amplitude
average or or
speed oflow
low frequency,
frequency,
the therate
the rateof
particle increases ofincrease
with increase ofofthe
increasing the average particle
average
amplitude particle
and fre-
speed
speed is lower.
is
quency. lower.
At low Moreover,
Moreover,
amplitude the
the amplitude
or amplitude has
low frequency, a greater
has a the effect
greater
rate on the
effect particle
on the
of increase thespeed
ofparticle compared
speed
average com-
particle
tospeed
paredthe to
frequency. However, there are distortion points in this process,
is lower. Moreover, the amplitude has a greater effect on the particle speedphe-
the frequency. However, there are distortion points in this and this
process, phenomenon
and this com-
occurs
nomenon moretheoften
pared to occurs at high
more
frequency.often frequencies andare
at high frequencies
However, there high amplitudes
and highpoints
distortion because
in thisthe
amplitudes dynamic
because
process, the behavior
anddynamic
this phe-
of the particle
behavior
nomenon the is
ofoccurs more
particle difficult
is more
more often to predict
difficult
at high to under
predict
frequencies these
andunder
highconditions
these because
conditions
amplitudes because ofthe
thedynamic
because strong
of the
nonlinearity
strong of
nonlinearity
behavior the dynamic system.
of theisdynamic
of the particle system.to predict under these conditions because of the
more difficult
strong nonlinearity of the dynamic system.
Separations 2022, 9, x FOR PEER REVIEW 11 of 15

Separations 2022, 9, x FOR PEER REVIEW 11 of 15

Separations 2022, 9, 216 11 of 14

Figure 7. Average speed with A = 1–10 mm, f = 8–26 Hz, θ = 10°, and a = 90°.

◦ and a = 90◦ .
Figure 7. Average
7.
FigureThe figures speed with A
with
and material
speed 1–10 mm, ffabove
A ==presented == 8–26
8–26discuss
Hz, θθ == 10°,
Hz, 10
the,and
influence
a = 90°.of common screening
parameters (amplitude and frequency) on the average particle speed affected, further-
The
The figures and and material presented above above discuss the influence of common screening
screening
more, byfigures
an inclination material
angle,presented
friction coefficient, discuss the influence
and vibration of common
direction angle. These pa-
parameters
parameters (amplitude
(amplitude and frequency)
and infrequency) on the average particle speed affected, furthermore,
rameters might not change the actualon the average
screening particle
process, speed
but they haveaffected,
a largefurther-
impact
by an inclination
more, angle, friction coefficient, and vibration direction angle. These parameters
on the average particle speed. Figure 8 shows the influence of the inclination angle on pa-
by an inclination angle, friction coefficient, and vibration direction angle. These the
might
rameters notmight
change notinchange
the actual screening process, but theybut have a large impact impact
on the
average particle speed in 5 sinattheA =actual
4 mm,screening
f = 16 Hz,process,
θ = 6–20°, and theya =have
90°. aInlarge
general, the
average
on particleparticle
the average speed. Figure 8Figure shows the influence of the inclination angle onanglethe average
average particle speed speed. increases with8 increasing
shows the influence
inclinationofangle.
the inclination
At the same on the
time, we
particle speed in 5 s at A = 4 mm, f = 16 Hz, θ = 6–20◦ , and a = 90◦ . In general, the average
average
found that particle
a suddenspeedincrease
in 5 s atoccurs
A = 4 mm,at θ f= =14° 16 because
Hz, θ = 6–20°,
particle and a = 90°.
motion onInthegeneral,
screening the
particle speed increases with increasing inclination angle. At the same time, we found
average
surface no particle
longer speed increases
converges, andwith14° increasing
is the convergenceinclination angle. At
threshold in the
thissame time,Fur-
instance. we
that a sudden increase occurs at θ = 14◦ because particle motion on the screening surface
found that a sudden increase occurs at θ = 14° because
thermore, we compare the convergence threshold of particle speed for different ampli- particle motion on the screening
no longer converges, and 14◦ is the convergence threshold in this instance. Furthermore,
surface
tudes and no frequencies.
longer converges, Figureand 14° isthe
9 shows theparticle
convergencespeed in threshold
50 s at Ain= 5this mm, instance.
θ = 14°,Fur- 15°,
we compare the convergence threshold of particle speed for different amplitudes and
thermore,
a = 90°, and wef =compare
16, 24, 32 theHz.convergence
The speedthresholdof the particle of particle
on thespeed for
screening ◦ different
surface
◦ ampli-
always
◦
frequencies. Figure 9 shows the particle speed in 50 s at A = 5 mm, θ = 14 , 15 , a = 90 , and
tudes and frequencies. Figure 9 shows the particle speed in particle
50 s at Ano = 5 longer
mm, θ converges
= 14°, 15°,
fconverges
= 16, 24, 32with Hz. increasing
The speed of frequency.
the particle The onmotion of the
the screening surface always converges with
awhen
= 90°,the and f = 16, 24,angle
inclination 32 Hz. The speed
increases by of These
1°. the particle
results onshow
the screening
that the surface
speed of always
conver-
increasing frequency. The motion of the particle no longer converges when the inclination
converges
genceincreases with
of particles increasing frequency. The motion of the particle no longer converges
angle by is
1◦ .not strongly
These results correlated
show thatwith frequency,
the speed although an
of convergence of increase
particles in fre-
is not
when
quency the inclination
leads to an angle in
increase increases
speed. bycontrast,
In 1°. These theresults showangle
inclination that thehas speed
a great ofinfluence
conver-
strongly correlated with frequency, although an increase in frequency leads to an increase in
gence
on theof
speed. In particles
contrast, is
convergence thenot strongly
threshold
inclination correlated
of the
angle has a with
particle speed.
great frequency,
Figure 10
influence onalthough
shows theanparticle
the convergence increase in fre-
speed
threshold in
quency
50 s at A leads
= 4, to
7, an
10 increase
mm, θ = in
14, speed.
15°, a In
= contrast,
90°, and f =the
16 inclination
Hz. This
of the particle speed. Figure 10 shows the particle speed in 50 s at A = 4, 7, 10 mm, θ = 14, angle
figure has a great
illustrates influence
that ampli-
on
15 ◦ the
tude =convergence
, a does90◦affect
, and fthe = threshold
particle
16 Hz. This of the
speed, particle
but
figure speed.that
the amplitude
illustrates Figure
and10 shows
speed
amplitude theaffect
particle
convergence
does speed
thethreshold
particle in
50 s at A
are slightly
speed, =
but the4, 7, 10 mm,
correlated. θ = 14,
amplitudeSimilarly, 15°,
and speed a = 90°,
weconvergence and
found that the f = 16 Hz. This
vibration
threshold figure
are directionillustrates that ampli-
angle has Similarly,
slightly correlated. little cor-
tude
we does
relation
found withaffect
that thethe
the particle direction
convergence
vibration speed,
of thebut the amplitude
particle
angle speed.
has and speedwith
little correlation convergence threshold
the convergence of
are slightly
the particle speed. correlated. Similarly, we found that the vibration direction angle has little cor-
relation with the convergence of the particle speed.

Figure 8.
Figure 8. Average
Averageparticle
particlespeed
speedin
in55ssat
atAA== 44 mm,
mm, ff = 16 Hz, θ 6–20◦ , and
θ = 6–20°, 90◦ .
and a == 90°.

Figure 8. Average particle speed in 5 s at A = 4 mm, f = 16 Hz, θ = 6–20°, and a = 90°.

Separations 2022, 9, x FOR PEER REVIEW 12 of 15

Separations 2022, 9, x FOR PEER REVIEW 12 of 15

Separations 2022, 9, x FOR PEER REVIEW 12 of 15
Separations 2022, 9, 216 12 of 14

Figure 9. Particle speed in 50 s at A = 5 mm, θ = 14, 15°, a = 90°, and f = 16, 24, 32 Hz.
Figure 9. Particle speed in 50 s at A = 5 mm, θ = 14, 15°,◦ a = 90°,◦and f = 16, 24, 32 Hz.
Figure Particlespeed
Figure9.9.Particle speedin 50ssatatAA==55mm,
in50 mm,θθ==14, 15 ,aa==90°,
14,15°, , andf =f 16,
90 and = 16,
24,24,
3232 Hz.
Hz.

Figure 10. Particle speed in 50 s at A = 4, 7, 10 mm, θ = 14, 15°, a = 90°, and f = 16 Hz.
Figure
Figure10.
Figure 10.Particle
10. Particlespeed
Particle speedin
speed in50
in 50sssat
50 atat
AA==4,4,7,7,10
A = 4, 7, 10 10mm,
mm,θ
mm,
= 14,
θθ==the
14, 15◦ ,aaa===90°,
14,15°,
15°, 90◦and
90°, , andff ==f 16
and
Hz.
= 16
16 Hz.
Hz.
Furthermore, the friction coefficient affects speed convergence threshold of the
particle. Figure 11 shows
Furthermore,
Furthermore, the the convergence
thefriction
friction coefficient condition
coefficientaffects
affectsthe theof the inclination
speed
speed convergence
convergence angle and friction
threshold
threshold ofthe
of the
Furthermore,
coefficient at A = 5 the friction
mm, a = 90°, coefficient
and f = 16 affects
Hz. The the speed
particle convergence
speed converges threshold
when of the
thefriction
value
particle.Figure
particle. Figure11 11shows
showsthe theconvergence
convergencecondition conditionof ofthe
theinclination
inclinationangle
angleandand friction
particle.
of Figure 11angle
the inclination showsand the convergence condition of the
the curve,
inclination angle and friction
coefficient
coefficient atatAA= = 5 mm,
5 mm, =friction
a =a90°, 90and fcoefficient
◦ , and
= 16f = 16The
Hz. Hz.is particle
below
The particle
speedspeed andconverges
convergesthe particle
when when
the speedthe
value
coefficient
diverges to at A = 5 mm,
infinity when a =the
90°,value
and fis= 16 Hz. the
above Thecurve.particle speed converges
Meanwhile, the when
curve of the
the value
con-
value of the inclination angle and friction coefficient is below the
of the inclination angle and friction coefficient is below the curve, and the particle speed curve, and the particle
of the inclination
vergence condition angle and friction coefficient isthebelow the curve, and the of particle speed
speed diverges
diverges to is
to infinity approximately
infinity
when thewhen
value theislinear,
aboveand
value istheabove convergence
theMeanwhile,
curve. threshold
the curvethe
curve. Meanwhile, the inclina-
curve
of the of
con-
diverges
tion
the angle to infinity
increases
convergence when
with
condition the
the value
friction
is is above
coefficient.
approximately the
linear,curve.
and Meanwhile,
the the
convergence curve of
thresholdthe con-
of the
vergence condition is approximately linear, and the convergence threshold of the inclina-
vergence
inclinationcondition is approximately
angle increases linear, and the convergence threshold of the inclina-
tion angle increases with thewith the friction
friction coefficient.
coefficient.
tion angle increases with the friction coefficient.

Figure 11. Convergence curve of the inclination angle and friction coefficient at A = 5 mm, a = 90◦ ,
and f = 16 Hz.
Separations 2022, 9, x FOR PEER REVIEW 13 of 15

Figure 11. Convergence curve of the inclination angle and friction coefficient at A = 5 mm, a = 90°,
Separations 2022, 9, 216 and f = 16 Hz. 13 of 14

Finally, the effect of the vibration direction angle on the average particle speed is
discussed.
Finally,Figure 12 shows
the effect of thethevibration
relation between
directionthe vibration
angle on the direction
averageangle andspeed
particle averageis
particle speed at A = 4 mm, f = 16 Hz, θ = 10°, a = 90°, and fc = 0.5. The average speed is
discussed. Figure 12 shows the relation between the vibration direction angle and average
minimalspeed
particle (Va = 0.3637
at A = m/s)
4 mm, when
f = 16theHz,
vibration
θ = 10◦ ,direction angle
a = 90◦ , and fc =is0.5.
90°.TheWhen the vibration
average speed is
direction angle starts to tilt with the vertical direction of the
minimal (Va = 0.3637 m/s) when the vibration direction angle is 90 . When thescreen ◦ surface, thevibration
average
speed increases significantly and reaches the maximum speed (V a = 0.7276 m/s) at a = 65°.
direction angle starts to tilt with the vertical direction of the screen surface, the average
The maximum
speed increasesaverage speedand
significantly is almost
reachestwice the minimum
the maximum speedaverage speed,m/s)
(Va = 0.7276 whichat aproves
= 65◦ .
that maximum
The the vibration direction
average speedangle has a large
is almost twiceinfluence
the minimumon theaverage
averagespeed,
speed.which
In theproves
actual
screening process, the granular material flow speed can be changed by adjusting
that the vibration direction angle has a large influence on the average speed. In the actual the vi-
bration direction angle. The average velocity begins to decrease and
screening process, the granular material flow speed can be changed by adjusting the oscillates with a fur-
ther reduction
vibration of theangle.
direction vibration
The direction angle. begins to decrease and oscillates with a
average velocity
further reduction of the vibration direction angle.

Figure 12. Relation between

between the the vibration
vibration direction
directionangle
angleand
andparticle
particleaverage
averagespeed
speedatatAA==4 4mm,
mm,f
f= =1616Hz,
Hz,θ θ= =10°,
10◦and fc =fc0.5.
, and = 0.5.

4. Conclusions
4. Conclusions
In this research,
In this research, the
thedynamic
dynamicequations
equationsofofa particle
a particleonon a screening
a screening surface
surface areare es-
estab-
tablished based on the elastoplastic contact model of round particles, and the
lished based on the elastoplastic contact model of round particles, and the dynamic be- dynamic
behavior is discussed.
havior is discussed. WeWe provide
provide a theoretical
a theoretical method
method forfor studying
studying thethe dynamic
dynamic behav-
behavior
ior
of coal granular material on a vibrating screen. This method can be applied to other other
of coal granular material on a vibrating screen. This method can be applied to gran-
granular materials
ular materials as well.
as well. The main
The main conclusions
conclusions are
are as as follows:
follows:
(1) The P–VS
(1) The P–VS system
system is
is strongly
strongly nonlinear.
nonlinear. A A small
small change
change inin parameters, such as
parameters, such as the
the
initial falling height and radius of the particle, will significantly affect the trajectory
initial falling height and radius of the particle, will significantly affect the trajectory of
the particle.
of the particle.
(2) In the
(2) In the normal
normal direction
direction of
of the
the vibrating
vibrating screen,
screen, thethe P–VS
P–VS motion
motion is
is quasiperiodic
quasiperiodic atat
low frequencies. With increasing frequency or amplitude,
low frequencies. With increasing frequency or amplitude, the motion of the motion of the
the P–VS
P–VS
system becomes chaotic, and its Poincaré map becomes petal-shaped. In addition, the
system becomes chaotic, and its Poincaré map becomes petal-shaped. In addition, the
number of petals increases at the mutation of the bifurcation diagram.
number of petals increases at the mutation of the bifurcation diagram.
(3) An increase in frequency, amplitude and inclination angle and a decrease in friction
(3) An increase in frequency, amplitude and inclination angle and a decrease in friction
coefficient lead to an increase in particle speed along the screen surface. In addition,
coefficient lead to an increase in particle speed along the screen surface. In addition,
the particle speed reaches a maximum when the vibration direction angle is 65◦ .
the particle speed reaches a maximum when the vibration direction angle is 65°.
(4) The divergence and convergence of particle motion along the screening surface are only
(4) The divergence and convergence of particle motion along the screening surface are
affected by the inclination angle and friction coefficient for the granular coal material.
only affected by the inclination angle and friction coefficient for the granular coal
In the future, the particle motion under optimal screening efficiency with different
material.
granular materials can be studied with our model. Whether the trajectories of different
In the future, the particle motion under optimal screening efficiency with different
granular materials under optimal screening are similar merits verification.
granular materials can be studied with our model. Whether the trajectories of different
granular
Author materials under
Contributions: optimal screening
Conceptualization, areC.L.;
D.H. and similar merits verification.
methodology, D.H. and C.L.; software, D.H.;
validation, D.H. and S.L.; formal analysis, D.H. and S.L.; investigation, D.H.; resources, D.H. and C.L.;
data curation, D.H. and S.L.; writing—original draft preparation, D.H. and S.L.; writing—review and
editing, D.H. and C.L.; visualization, D.H.; supervision, D.H.; project administration, C.L.; funding
acquisition, C.L. All authors have read and agreed to the published version of the manuscript.
Separations 2022, 9, 216 14 of 14

Funding: The research presented in this work was funded by the National Natural Science Founda-
tion of China [51775544] and Xuzhou XCMG Mining Machinery Co, Ltd. [2020050049-07].
Data Availability Statement: The data is availabile to contact the author ([email protected]).
Conflicts of Interest: The authors declare no conflict of interest.

References
1. Peng, L.; Wang, Z.; Ma, W.; Chen, X.; Zhao, Y.; Liu, C. Dynamic influence of screening coals on a vibrating screen. Fuel 2018, 216, 484–493.
[CrossRef]
2. Chen, Z.; Li, Z.; Xia, H.; Tong, X. Performance optimization of the elliptically vibrating screen with a hybrid MACO-GBDT
algorithm. Particuology 2021, 56, 193–206. [CrossRef]
3. Jiang, H.; Qiao, J.; Zhao, Y.; Duan, C.; Luo, Z.; Liu, C.; Yang, Y.; He, J.; Zhao, L.; Pan, M. Evolution process and regulation of
particle kinematics and spatial distribution driven by exciting parameters during variable-amplitude screening. Powder Technol.
2018, 330, 292–303. [CrossRef]
4. Horenstein, M.N.; Mazumder, M.; Sumner, R.C. Predicting particle trajectories on an electrodynamic screen—Theory and
experiment. J. Electrost. 2013, 71, 185–188. [CrossRef]
5. Jiang, Y.-Z.; He, K.-F.; Dong, Y.-L.; Yang, D.-l.; Sun, W. Influence of Load Weight on Dynamic Response of Vibrating Screen. Shock
Vib. 2019, 2019, 4232730. [CrossRef]
6. Zhao, L.; Zhao, Y.; Bao, C.; Hou, Q.; Yu, A. Optimisation of a circularly vibrating screen based on DEM simulation and Taguchi
orthogonal experimental design. Powder Technol. 2017, 310, 307–317. [CrossRef]
7. Pan, J.W.; Li, J.; Hong, G.Y.; Bai, J. A mapping discrete element method for nonlinear dynamics of vibrating plate-particle coupling
system. Powder Technol. 2021, 385, 478–489. (In English) [CrossRef]
8. Jiang, H.; Yu, S.; Huang, L.; Zhao, Y.; Cao, X. Kinematics and mechanism of rigid-flex elastic screening for moist coal under
disequilibrium excitation. Int. J. Coal Prep. Util. 2020, 42, 1724–1739. [CrossRef]
9. Wang, L.; Ding, Z.; Meng, S.; Zhao, H.; Song, H. Kinematics and dynamics of a particle on a non-simple harmonic vibrating
screen. Particuology 2017, 32, 167–177. [CrossRef]
10. Yang, Y.; Wan, L.R. Study on the Vibroimpact Response of the Particle Elastic Impact on the Metal Plate. Shock Vib. 2019, 2019,
1–13. (In English) [CrossRef]
11. Dong, J.; Fang, J.; Pan, J.; Hong, G.; Li, J. Dynamic model of vibrating plate coupled with a granule bed. Chaos Solitons Fractals
2022, 156, 111857. [CrossRef]
12. GFF, G. Hertz’s Miscellaneous Papers. Nature 1896, 55, 6–9.
13. Hunt, K.; Crossley, E. Coefficient of restitution interpreted as damping in vibroimpact. J. Appl. Mech. 1975, 42, 440–445. [CrossRef]
14. Thornton, C. Coefficient of Restitution for Collinear Collisions of Elastic-Perfectly Plastic Spheres. J. Appl. Mech. 1997, 64, 383–386.
[CrossRef]
15. Jian, B.; Hu, G.M.; Fang, Z.Q.; Zhou, H.J.; Xia, R. A normal contact force approach for viscoelastic spheres of the same material.
Powder Technol. 2019, 350, 51–61. [CrossRef]
16. Safaeifar, H.; Farshidianfar, A. A new model of the contact force for the collision between two solid bodies. Multibody Syst. Dyn.
2020, 50, 233–257. [CrossRef]
17. Shen, Y.; Xiang, D.; Wang, X.; Jiang, L.; Wei, Y. A contact force model considering constant external forces for impact analysis in
multibody dynamics. Multibody Syst. Dyn. 2018, 44, 397–419. [CrossRef]
18. Vu-Quoc, L.; Zhang, X.; Lesburg, L. Normal and tangential force–displacement relations for frictional elasto-plastic contact of
spheres. Int. J. Solids Struct. 2001, 38, 6455–6489. [CrossRef]
19. Yang, Y.; Cheng, Y.M. A fractal model of contact force distribution and the unified coordination distribution for crushable granular
materials under confined compression. Powder Technol. 2015, 279, 1–9. (In English) [CrossRef]
20. Ye, Y.; Zeng, Y.W.; Chen, X.; Sun, H.Q.; Ma, W.J.; Peng, Z.X. Development of a viscoelastoplastic contact model for the size- and
velocity-dependent normal restitution coefficient of a rock sphere upon impact. Comput. Geotech. 2021, 132, 104014. (In English)
[CrossRef]
21. Yu, J.; Chu, J.; Li, Y.; Guan, L. An improved compliant contact force model using a piecewise function for impact analysis in
multibody dynamics. Proc. Inst. Mech. Eng. Part K J. Multi-Body Dyn. 2020, 234, 424–432. [CrossRef]
22. Zhang, J.; Li, W.; Zhao, L.; He, G. A continuous contact force model for impact analysis in multibody dynamics. Mech. Mach.
Theory 2020, 153, 103946. [CrossRef]
23. Zhao, P.; Liu, J.; Li, Y.; Wu, C. A spring-damping contact force model considering normal friction for impact analysis. Nonlinear
Dyn. 2021, 105, 1437–1457. [CrossRef]
24. da Silva, M.R.; Marques, F.; da Silva, M.T.; Flores, P. A compendium of contact force models inspired by Hunt and Crossley’s
cornerstone work. Mech. Mach. Theory 2022, 167, 104501. (In English) [CrossRef]
25. Ji, S.; Liu, L. Contact Force Models for Granular Materials. In Computational Granular Mechanics and Its Engineering Applications;
Ji, S., Liu, L., Eds.; Springer: Singapore, 2020; pp. 51–96.
Reproduced with permission of copyright owner. Further reproduction
prohibited without permission.