Extension of the bouncing ball model to a vibratory conveying system

Various special effects occur during the operation of vibratory conveyors, e.g., multiple feeding velocities at the same excitation amplitude or so-called microthrows. In this work, a model for the simulation and prediction of the behavior of such a conveying system is presented. The simulation model is based on the bouncing ball model which is known from literature. The introduced impact law is coupled in horizontal direction by a frictional force which enables modeling a feeding process. The mentioned effect of multiple feeding velocities is studied with the developed simulation model. For the estimation of the critical excitation amplitude where a second feeding velocity appears, an analytical approach is developed. The corresponding feeding velocity can also be calculated with this approach. Moreover, the sensitivity of the initial conditions is investigated and criteria for the estimation are found. These can be applied to optimally adjust the conveyor in practice. Furthermore, the effects of microthrows are studied and analytical formulas for the estimation of characteristic values of the microthrows are derived. The dragging process following a sequence of microthrows is also investigated. All the developed formulas are validated by the simulation model.


Introduction
The demands on the automation industry are becoming higher due to increasing quality and customer requirements, which is why many companies are presented with enormous challenges.Therefore, increasing the efficiency of the respective systems is also becoming increasingly important.However, this requires a deep understanding of the system in order to be able to make improvements.Especially in the field of conveyor technology, there are many unanswered questions concerning the nonlinear behavior of the conveyor.Some of these questions are addressed in this paper to enable a deeper understanding of the system.
There are already many approaches for the simulation of vibratory conveying systems.However, there appear many special effects which have to be modeled in a simulation tool.In [1], Gravenkötter investigates new phenomena in vibratory conveyors by examining adhesive forces and comparing measurement to improved model results.The work of Nendel [2] introduces methods for analyzing two-dimensional movement patterns of vibratory conveyors and examining the use of vibration isolators to reduce dynamism transfer.In [3], a numerical 2D model is developed to investigate the dynamic behavior of a feeding part in a vibratory bowl feeder.It is found that the mean conveying velocity depends on the part shape, vibration angle and coefficient of friction.However, it can exhibit both periodic and chaotic behavior.The dissertation thesis of Hofmann [4] aims to develop a simulationbased approach for designing sorting chutes in vibratory conveyors.The goal is to accurately model the essential physical effects of conveying and sorting processes.The simulation results will be used to enhance the practical development process of conveyors.In the study of Kobayakawa et al. [5], particle saltation on an obliquely oscillating plate using a mass-point model is simulated.It is founded that particles with high restitution bounce forward and backward repeatedly, while particles with low restitution only bounce forward due to the phase angle of the oscillating plate and impulse during particle collision.In the work of Bednarski [6], the suitability of existing vibratory transport models for the metallurgical industry is examined.The analysis reveals that these models neglect the influence of air flow through the feed layer, leading to reduced transport velocity when handling dusty and powdery fractions.To address this issue, a mathematical and simulation model is formulated to incorporate the air flow effect.Experimental data is used to validate the model.In [7], the performance of a vibratory conveyor is optimized.Furthermore, the vibration transmission is reduced through dynamic model development, comprehensive analysis and optimization.In [8], a method for modeling and simulating linear feeding systems in industrial automation is presented.Czubak's study [9] focuses on analyzing new solutions for vibratory conveyors, specifically for transporting loose materials or small objects at variable velocities.The investigation explores ways to interrupt material flow without shutting down the drive.Analytical and simulation analyses of the conveyor's performance show promising results, providing a foundation for further laboratory investigations.In [10], the motion of a small rigid particle on an oscillating inclined conveying tray is simulated.The longitudinal travel due to differential friction with small oscillation amplitudes is observed, while larger amplitudes lead to periodic sliding phases.Chaotic bounc-ing behavior is observed at even larger amplitudes, but periodic bouncing can be achieved when the particle is dropped from a specific height with certain parameters.In the research conducted by Surowka [11], a novel vibratory conveyor for precise material dosage is studied.The study analyzes the transport capabilities of this conveyor in the resonance zone using analytical and simulation methods.An optimal operating point for the system is identified, and simulation results determine transport velocities at different excitation frequencies.The findings are confirmed using a specially designed industrial conveyor, based on the analytical and simulation investigations.
A main challenge in modeling a conveying process is the contact between part and conveyor.In [12], various methodologies are used to develop impact models for colliding bodies.A comparison between different models is also made.A further comparison between several contact models is performed in [13].Additionally, nonlinearities and micro-slip effects are treated therein.The paper of Michalczyk [14] confirms the constancy of the coefficient of restitution in collisions involving material damping.However, he suggests that the coefficient of restitution may depend on the density of energy flux in collisions of different nature, including eccentric collisions.He proposes the concept of describing the coefficient of restitution as a random function of the density of the energy flux.The work of Hertz from 1881 presents fundamental considerations on the contact between two bodies [15].Further basic analyses on contact mechanics are presented in [16,17].In [18], a simplified model of a conveyor is developed which allows to study the movement of a point mass on a moving plate with numerical experiments.The formulation of the contact model is based on the works in [15][16][17].
If only the vertical movement of material in a conveying process is considered, this corresponds to the bouncing ball model where a point mass bounces on a moving plate and switches between contact and free fall.A lot of literature can be found on this.The work of Guckenheimer from 1984 [19] applies techniques from dynamical systems and bifurcation theory to the study of nonlinear oscillations, with an emphasis on the geometrical and topological properties of solutions.It is a fundamental work for the understanding of nonlinear theory.In [20], a horseshoe map of the bouncing ball problem is created for the evidence of chaotic behavior.In [21], Tufillaro focuses the effect of period doubling of the bouncing ball model.Furthermore, the influence of frequency and amplitude in order to identify the separation of periodic orbits in parameter space is investigated.In [22], the mechanical properties of a ball that is dropped vertically are related to the coefficient of restitution.In many areas of application, the stability of rigid bodies plays an essential role, for example, in satellite technology [23].Another, more theoretical, field of application is the analysis of the stability of fixed points [24].There are many works which investigate the solution of fixed point problems in general, e.g., [25,26].The stability analysis also finds application in the bouncing ball problem [19].A complete derivation of the stability criterion for the bouncing ball problem is presented in [27,28].
In this work, a simulation model is presented which allows to predict the behavior of a conveying system.The model is based on the bouncing ball formulation which describes the behavior of a point mass on a moving plate in vertical direction.The parameter which controls the impact is the coefficient of restitution.It is defined as the ratio of the normal velocity after vs. before the impact.Therefore, it is a measure for the energy dissipation during the impact.The scientific novelty of this manuscript is the coupling of the wellknown motion in vertical and horizontal direction.This coupling is performed by a frictional force and enables to simulate the complex feeding process, see Sect. 2. Furthermore, the simulation model allows to model the nonlinear effects, occurring during a feeding process.Also the chaotic behavior of the system can be predicted and analyzed with the developed simulation model.
A main effect observed in measurements of conveying systems is the occurrence of multiple feeding velocities at the same excitation amplitude.This effect is very sensitive and there is still a lack of methods to predict this behavior.In Sect.3, an analytical method is developed which allows to estimate the critical amplitude.This method is again based on the bouncing ball formulation and is extended to feeding direction.Furthermore, this enables to transfer the results to a conveying process and yields an estimation for a maximum conveying velocity.The analytical formula allows to evaluate this for several parameter settings and gives a generalization of the complex process.However, these results can be applied to practical application which enables a massive improvement of the performance of the conveyor.A further main effect is the sensitivity of the second state of motion which appears above the predicted critical amplitude.For the analysis of the sen-sitivity, a stability criterion of the bouncing ball is used.This method is also presented in Sect.3.
A further main effect in vibratory conveying systems is the so-called microthrows.These are characterized by the fact that the relative speed is very low and the part is dragged along if the relative speed is converged to zero.In Sect.4, general approaches and formulas to understand these complex motions are derived which is a further scientific novelty.In [29], Vogel and Linz say that there is no analytical expression until this time.Furthermore, the bouncing ball model is extended to describe the dragging process after the microthrows which it is not possible with the original equations.However, this extension is done for both, the vertical and the horizontal direction.
The application of the presented methods to the conveying system is discussed in Sect. 5.The different states of motion including microthrows are discussed and the developed analytical approaches are shown on an example.Furthermore, the critical excitation amplitude where the second state of motion appears is investigated with regard to dependencies and sensitivity.In addition, remarks about the transferability into practice are given throughout the different sections.

Contact model and the extended bouncing ball problem
In this section, a model is developed which describes a conveying process of a point mass on a periodically moving conveyor.Especially the impact in vertical direction and the coupling between vertical and horizontal direction are discussed.The vertical direction is also called the normal or z-direction, the horizontal direction the tangential or x-direction.In Fig. 1, the described situation is shown.The variables v N and v T describe the velocity of the part in normal and tangential direction.The velocities of the conveyor are described with w N and w T .The notation v denotes the velocity before the impact and v denotes the velocity after the impact.However, it is assumed that the velocity of the conveyor does not change during the impact.
The derivation of the contact model is based on Newton's law in normal and tangential direction with (2)

Impact law for normal direction
For modeling the instantaneous contact situation, an impact law in normal direction is introduced, see Fig. 2.
It is assumed that the duration of the impact is equal to ε which is close to zero.Furthermore, the norm of the force is chosen with S ε , where S denotes the momentum of the impact.
The consequence of this normalization is that the integral of the force in the normal direction over time during contact is equal to the momentum S, see Eq. (3).
Therefore, the integration of Eq. ( 1) results in which is equivalent to In Eq. ( 5), a connection between the impact S and the velocity difference before and after the impact is found.During the contact from t = 0 to t = ε, energy is dissipated because of the damping of the contact.This damping is described with the coefficient of restitution c r which describes the ratio of the relative velocity after and before the impact [12].With the assumption that the velocity of the conveyor does not change during the contact, the velocity after the impact v N can be described with Inserting Eq. ( 6) in Eq. ( 5) results in In Eq. ( 7), the unknown velocity after the impact is described with the coefficient of restitution c r .

Impact law for tangential direction
In order to model the coupling of the vertical component with the tangential component, a coefficient of sliding friction μ is introduced.Furthermore, the velocity difference of conveyor and part in feeding direction has to be considered by computing the tangential force f T with Inserting Eq. ( 8) in Eq. ( 2) results in and the differential equation in Eq. ( 9) can be solved.The term (w T − v T ) can be positive or negative which leads to a change of sign, see Eq. ( 9).Therefore, the tangential velocity of the part during the contact can be described with a time dependent function, see Eq. ( 12).
However, a few cases have to be distinguished which are discussed in the following.In Fig. 3, two possible cases (a) and (b) are shown which may occur if the velocity of the part before the impact v T is lower than the velocity of the conveyor w T (a) v T < w T and v T < w T : In case (a), the tangential velocity of the part before and after the impact is lower than the tangential velocity of the conveyor, which results in sign(w T − v T ) = 1.Therefore, the velocity of the part after the impact can be calculated with In case (b), the tangential velocity of the part before the impact is lower than the tangential velocity of the conveyor.After the impact, it would be greater than the tangential velocity of the conveyor which is physically not possible.Therefore, the velocity after the impact is equal to the velocity of the conveyor.
In Fig. 4, two situations are illustrated which may occur if the velocity of the part before the impact v T is larger than the velocity of the conveyor w T .
(c) v T > w T and v T > w T : In case (c), the tangential velocities of the part before and after the impact are larger than the tangential velocity of the conveyor, which results in sign(w T − v T ) = −1.Therefore, the velocity of the part after the impact results in In case (d), the velocity of the part before the impact is larger than the velocity of the conveyor.Applying Eq. ( 15), it would be lower than w T after the impact which is physically not possible again, compare case (b).Therefore, the velocity after the impact is again equal to the velocity of the conveyor.v T = w T (16)

Poincaré map for modeling the conveying process
Next, we introduce a Poincaré map which describes the state of the conveying system immediately after impact k + 1 as a function of the state immediately after impact k.The state after impact k is given by the dimensionless impact time θ k = ωt k , where ω denotes the excitation frequency of the conveyor, and by the normal and tangential velocities v N ,k and v T,k .So, we are looking for the map This map can be interpreted as an extension of the bouncing ball model, as, e.g., described in [27].In Fig. 5, the conveying situation with the extended model is shown.During the phases without contact with the oscillating plate, the transported mass moves on a flight parabola.On the other hand, each contact causes a sudden change of the particle velocity according to the impact laws in tangential and in normal direction as described above.
Assuming a sinusoidal motion of the conveyor plate with the frequency ω and the vertical displacement amplitude A, the impact time θ k+1 can be computed from the condition that the vertical position of the plate and the mass point must be equal at each impact time.The vertical conveyor position is given by z c (t) = A sin(ωt) and the vertical position of the transported particle after impact k reads From z c (t k+1 ) = z p (t k+1 ) an by substituting t k+1 = θ k+1 /ω, we obtain This relation can be solved for θ k+1 , but note that one must find the smallest positive root.In the next step, the normal velocity v N ,k+1 can be computed from Eq. ( 6) Finally, the tangential velocity v T,k+1 is determined by one of Eqs.[13][14][15][16], where v T = v T,k .Moreover, if x c (t) = B sin(ωt) denotes the plate motion in horizontal direction, we identify w T = ẋc (t k+1 ) = ωB cos(θ k+1 ).The momentum of the impact k + 1 introduced in Eq. (7), which is also needed for computing v T,k+1 , is given by Using this, the tangenial velocity v T,k+1 according Eqs. ( 13)-( 16) can be expressed in compact form as follows: 3 One-periodic motions In a vibratory conveying process, several states of motion can be distinguished.In one state of motion, the part is periodically hopping on the conveyor.This motion appears only above a critical excitation amplitude [18].In Sect.3.1, an approach based on the bouncing ball model is developed for the prediction of this amplitude and the feeding velocity related with this motion.Besides the prediction of the critical amplitude, the stability of the periodic motion is investigated in Sect.3.2.

Conditions for the appearance of one-periodic motions
For finding sufficient conditions for the appearance of periodically hopping motions, we look for a oneperiodic solution of Eqs.19 and 20.Mathematically speaking, the phase θ k should continue 2π -periodically with where n denotes the number of periods the conveyor moves during the flight of the part.Hence, where θ * = mod (θ k , 2π) is a proper phase angle.Moreover, the velocity in normal direction v N ,k = v * N Fig. 5 Conveying process of point mass m on moving plate has to be the same at every point of impact.Inserting Eq. ( 23) into Eq.( 19), results in From this, the normal velocity v * N for a one-periodic motion can be computed as Inserting this velocity for v N ,k+1 and v N ,k in Eq. ( 20), yields the following phase angle for the one-periodic motion: Equation 26 shows that a real number for θ * only exists if cos(θ Hence, a one-periodic motion can only occur if The feeding velocity associated with the periodically hopping motion can be determined from the impact laws described in Eq. (22).Depending on the initial velocity in tangential direction, the feeding part is either accelerated or decelerated after each impact until it moves with the same velocity as the conveyor plate.Hence, in a steady state, the part is transported with the tangential plate velocity at time θ * , i.e., with w T = ωB cos(θ * ).Inserting Eq. ( 26) for θ * yields Note, that v * T depends only on the restitution coefficient c r , on the excitation frequency ω, and on the ratio B/A of the horizontal and normal excitation amplitude, i.e., only on the direction in which the plate is moving, but not on the amplitude of the motion itself, see also Eq. (67).

Stability analysis of the one-periodic motion
Next, a stability criterion for the one-periodic solution of the bouncing ball problem is derived.A similar derivation can be found in [28].
For that purpose, a small perturbation of the periodic solution Eqs. 24 and 25 is introduced by setting ) After inserting in Eqs.19 and 20 and linearizing with respect to the perturbations δθ k and δv k , the linear mapping is obtained, where the conditions for v * N and cos θ * shown above have been used.The one-periodic solution is stable if the absolute value of both eigenvalues of the mapping matrix is smaller than one.The characteristic equation for the eigenvalues reads where The solution of the characteristic equation is given by the eigenvalues satisfy since the coefficient of restitution is always smaller than one.Hence, the periodic orbit is stable for z within the range described in Eq. (34).However, if z 2 > 4c 2 r , i.e., if z > 2c r or z < −2c r the eigenvalues are real numbers which, for stability, must satisfy the conditions where the first condition must be investigated for z > 0 and the second for z < 0. In this cases, the periodic orbit is stable if Combining the stable domains for z from Eqs. 34 and 37 finally yields the stability condition Inserting the definition of z from Eq. ( 32), the stability condition reduces to If sin θ * , which must be positive, is substituted by √ 1 − cos 2 θ * and Eq. ( 26) is used for expressing cos θ * in terms of the system parameters, the range of the excitation amplitude for stability reads where the lower limit is equivalent to the lower bound of the amplitude of the one-periodic motion in Eq. ( 27).

Microthrows
In the bouncing ball problem, Eqs. ( 19) and (20), which govern the motion in normal direction, may converge to By inserting 19 and 20 it can be seen that these relations are satisfied.In this case, the time between two impacts tends to zero whereas the normal velocity converges to the conveyor velocity.We call this motion a sequence of microthrows.

Transport velocity during a sequence of microthrows
For analyzing a sequence of microthrows, we may approximate the conveyor motion between the impacts k and k + 1 by a parabola: Hence, at t = t k+1 and using again θ = ωt, the conveyor position is given by where Δθ k := θ k+1 −θ k is the (normalized) time for the kth throw.From inserting in Eq. ( 19), we then obtain Solving for Δθ k yields where is the difference between the particle and the conveyor velocity in normal direction after the kth impact, which cannot become negative.Therefore, Δθ k is positive, i.e., θ k+1 > θ k , only if g − Aω 2 sin θ k > 0. Hence, a sequence of microthrows can only occur if the phase angle satisfies the condition where Γ is the so-called driving factor.This factor describes the acceleration amplitude of the conveyor related to the gravity and is used in many works for describing the excitation as dimensionless number, see, e.g., [27,28].
The conveyor velocity after the kth impact results from differentiation of Eq. ( 42) Using Eq. ( 18), the relative impact velocity with which the particle hits the plate at t = t k+1 is given by where Eq. ( 44) has been inserted.Hence, after the impact, the difference velocity between particle and plate is given by Note that Δv N ,k is reduced after each impact by the factor c r as for a resting plate.For k → ∞, Δv N ,k converges to zero.To compute the decrease of the flight time Δθ k , we derive from Eq. (44) that If Eq. ( 46) holds for k and k + 1, the factor on the right side is positive.Since c r ≤ 1, the ratio Δθ k+1 /Δθ k is less than one if This yields a convergent series where Δθ k also converges to zero.Expanding the factor on the right side of Eq. (49) reads

and, hence,
By using the summation formula for infinite series, the total time for a sequence of microthrows can then be estimated as Next, we want to compute the mean horizontal velocity of the conveyed particle.Recall that v T,k denotes the horizontal velocity during the kth throw.The mean horizontal velocity for a sequence of microthrows is given by For computing v T,k , we have to determine the momentum S k for a sequence of microthrows.By setting cos θ k+1 ≈ cos θ k −Δθ k sin θ k , we obtain from Eq. ( 21) where Δθ k has been substituted from Eq. (44).If the conveyor is moving faster than the conveyed part throughout the whole sequence of microthrows, the horizontal velocity is given by Eq. ( 22) as Using Eq. ( 48) and the summation formula for series, we obtain Hence, the tangent velocity after the sequence of microthrows is given by provided that v T,∞ < ωB cos(θ 1 + Δθ ).Moreover, after inserting Eq. (55) in Eq. ( 52), the mean transport velocity reads (57)

The bouncing ball model for θ > θ ∞
Obviously, the bouncing ball model cannot describe what happens after a sequence of microthrows for which the Poincaré map Eq. ( 17) converges to a fixed point θ ∞ .However, the system will of course continue to move in some way.Hence, the question arises how the bouncing ball model can be extended to compute states for θ > θ ∞ .After a sequence of microthrows, the vertical motion of the conveyed part and the conveyor will be equal.Therefore, The normal force F N acting from the conveyor on the part results from Newton's second law m z p = F N −mg and is, hence, given by The coupled normal motion ends when F N becomes negative, i.e., when Then, the Poincaré map can be initiated again with θ 0 and v N ,0 = Aω cos θ 0 .For the initial tangent velocity v T,0 after coupled normal motion, we use Newton's second law in tangential direction under consideration of Eq. ( 8) with where θ ∞ and v T,∞ and denote the phase angle and the tangential velocity at the beginning of the coupled normal motion, i.e., at the end of the preceeding sequence of microthrows.

Results and discussion
A problem in practical applications is that several conveying speeds may occur at the same excitation amplitude.This effect can be observed on a test rig where the transport time of a single part is measured between two light barriers several times.It is assumed that the initial conditions are responsible for this phenomenon [30].
For reproducing this effect with a simplified model, the initial height is varied at several excitation amplitudes.
The mean feeding velocity is evaluated at a steady state [18].However, it has to be mentioned that the representation of the results in this paper does not correspond to a classical bifurcation diagram as it is used in bouncing ball theory [19,21,22].This is due to the output variable used, including the constant excitation amplitudes and the variation of the initial conditions.Therefore, the effect of period doubling before the transition to chaotic behavior [27] may not be detected.However, for showing the appearance of multiple states of motion at the same excitation amplitude this representation is suitable.
The excitation of a conveying system is typically given with a frequency f , an acceleration amplitude â and an angle of throw ϕ which denotes the angle between horizontal and vertical component.Therefore, the vertical displacement amplitude can be calculated with A = â sin(ϕ) ω 2 and the horizontal displacement amplitude with B = â cos(ϕ) ω 2 .The default set of excitation parameters is defined with the frequency f = 50 Hz and the angle of throw ϕ = 12 • .Furthermore, it has to be mentioned that in most cases the conveyor moves one period during the flight of the part.Therefore, by default, it is specified that n = 1, see Eq. (23).The default set of contact parameters is defined with the coefficient of restitution c r = 0.6 and the coefficient of sliding friction μ = 0.15.The mass of the part is m p = 0.7 × 10 −3 kg and it starts at a height of z 0 = 0 mm with an initial velocity v 0 = 0 mm/s.The code for the simulation of the conveying process is implemented in the programming language Python.The details for modeling several states of motion are discussed in the following.

States of motion
The main result analyzed in this section is the mean feeding velocity of the part in tangential direction in a Fig. 6 Multiple feeding velocities by varying the initial position z 0 at every acceleration amplitude â steady state motion.For an illustration of the possible states of motion, the initial position z 0 of the point mass is slightly varied at different acceleration amplitudes â, see above.The results of this analysis with the full simulation model including microthrow effects are shown in Fig. 6.The results consist of three different states of motion.The first state of motion results in a lower mean feeding velocity and appears below â < 60 m/s 2 .The second state of motion results in a higher mean feeding velocity which appears at â = 42 m/s 2 and ends at â = 62 m/s 2 .It has to be mentioned that this is the maximum stable velocity of a one-periodic motion with n = 1 which is studied in detail afterward in this section.Between â = 42 m/s 2 and â = 60 m/s 2 , there exist two states of motion at the same acceleration amplitude.The different feeding velocities are only due to the small variation of the initial condition z 0 .This phenomenon is studied in [18] and is further investigated in Sect.5.3.Above â = 62 m/s 2 , there is a chaotic state of motion [18].Next, the time signals of the respective state of motion are discussed.In Fig. 7, the time signals of the second state of motion at â = 50 m/s 2 can be seen.It is obvious that the feeding velocity in tangential direction is constant in the steady state motion where the point mass moves periodically in a hopping mode on the conveyor.The periodic signals can be generated with the Poincaré map in Eq. ( 17) without any extension.
In Sect.3, a method for predicting the feeding velocity of a one-periodic motion has been derived.The resulting formula is shown in Eq. (28) and yields a predicted feeding velocity of v T = 11.5 cm/s using the default set of parameters with c r = 0.6.The upper  7 shows that the predicted feeding velocity (orange) coincides with the simulation results (blue).A more detailed analysis of predicting the feeding velocity and the sensitivity of the initial conditions follows in Sect.5.2.
In Sect.4.1, it is shown that the relative normal velocity may approach zero which is called sequence of microthrows, see Eq. ( 48).It should be noted that these microthrows can still be simulated with the bouncing ball equations in Eqs.19 and 20.Only if the time of impact θ k has converged to θ ∞ , i.e., the relative normal velocity is zero, no further solution can be calculated.
Therefore, the simulation model is extended with the complements in Sect.4.2.It is shown that the velocity of the part is the same as the velocity of the conveyor after the sequence of microthrows.This is the case until the criterion Eq. ( 60) is fulfilled.If the part lifts off after this coupled normal motion, the bouncing ball equations in Eqs.19 and 20 can be used again.
In Fig. 8 the time signals of the first state of motion at â = 49 m/s 2 with the full simulation model are shown.The lower diagram shows the z-position of the part and the conveyor.It seems that the part does not lift off above t ≈ 0.1 s.However, it should be mentioned that the part lifts off slightly, which will be investigated in detail later.It is also noticeable that the conveying speed decreases even though the z-position of the part and the conveyor coincide.This can be explained by the fact that the transient process of the feeding velocity takes a longer time because it is coupled via friction force, see Eqs. (62) and ( 9).
In the following, the state of motion in Fig. 8 is analyzed step by step.As mentioned above, it starts with a sequence of microthrows.In Fig. 9, the time signals of the first sequence of microthrows can be seen, where the green lines denote the times of impact.In the lower graph it is shown that the time between the impacts decreases continuously and above t = 0.077 s, the impacts end.This is because the limit of Eqs.19 and 20 is reached, since θ k = θ ∞ .It can also be seen in the upper graph in Fig. 9 because there is no feeding velocity above t = 0.077 s.
The microthrow in Fig. 9 starts at t = 0.07518 s.The developed formula in Eq. (51) allows to predict the total time of the sequence of microthrows which yields Δt = 0.00186 s.This is a deviation of 6.79% compared to the simulation results in Fig. 9. Furthermore, the tangential velocity after the sequence at microthrows (t = 0.077 s) can be estimated with Eq. (56).This yields a feeding velocity of v T = 5.59 cm/s.The simulated feeding velocity at the end of the sequence of Fig. 10 Relative normal velocity during sequence of microthrows microthrows in Fig. 9 is v T,sim = 5.26 cm/s which represents a deviation of 6.34%.The mean feeding velocity during the sequence of microthrows can be predicted with Eq. ( 57) and results in v = 5.43 cm/s which differs by 4.69% from the simulated value, see Fig. 9. Summarizing, the formulas from Sect.4.2 considering the microthrows are validated by an example.The deviations results from the linearized assumptions in Eqs.51, 56 and 57.
In Fig. 10, the relative normal velocity can be seen.It is obvious that it is equal to zero at t = 0.077 s, so the bouncing ball model is not able to model a further impact, see Eq. (44).However, the system will continue to move in some way above the first sequence of microthrows.Therefore, the model is extended with the considerations in Sect.4.2.We know from Eq. (58) that the motion of the part after a sequence of microthrows coincides with the motion of the conveyor, which is called coupled normal motion.This coupled normal motion after the first sequence of microthrows for θ > θ ∞ is shown in Fig. 11, using the extended simulation model.In the domain where part and conveyor are in permanent contact (compare Eq. 58), only the z-position of the conveyor (orange) can be seen because it coincides with the z-position of the part.
The upper graph in Fig. 11 shows the feeding velocity of the part which is calculated with Eq. (62).As mentioned before, the transient process of the feeding velocity takes a longer time because it is coupled via friction force.However, between the times of impact at t > 0.08, the feeding velocity does not change with indicates a flight phase.
In Eq. (60), a criterion is derived which has to be fulfilled for the existence of a coupled normal motion.This means that the part can only lift off again if the acceleration of the conveyor is greater than the gravitational constant.In Fig. 12 it can be seen that this is the case at t = 0.083 s, which validates the simulation   11.Furthermore, the assumption that the part has to be in the flight phase between the impact lines in Fig. 11, is validated.However, this sequence of processes is repeated until the end of the simulation time, compare also Figs. 8 and 11.
In Fig. 13, the time signals of a chaotic state of motion at â = 80 m/s 2 are shown.It is obvious that there is no steady state and the system switches between movements in the hopping mode, sequences of microthrows and coupled normal motions.
A detailed discussion of the states of motion, especially the chaotic behavior, can be found in [18].

Prediction of multiple feeding velocities
In Fig. 14, the excitation amplitude is varied again with multiple initial conditions, compare Fig. 6.Furthermore, the maximum velocity of the conveyor in x-direction is shown.The described analysis is shown for two coefficients of restitution, c r = 0.4 (left) and c r = 0.7 (right).It can be seen that the appearance of the second solution differs due to the change in the coefficient of restitution.
In Fig. 15, the evaluation of Eq. ( 26) is shown for the two analyses with different coefficients of restitution.The critical excitation amplitude is described with cos(θ * ) and it is shown that it is equal to one at this acceleration amplitude where the second state of motion appears, compared with Fig. 14.The second state of motion does not occur if the excitation amplitude is less than the critical amplitude, see Eq. (27).
The critical amplitude where the second state of motion may appear, is solved generally in Eq. ( 27).However, for the application to the conveyor system, the solution in Eq. ( 27) has to be converted to an equivalent acceleration amplitude with where ϕ denotes the angle of throw between horizontal and vertical component.Therefore, the critical acceleration amplitude â reads â = πgn sin(ϕ) From Eq. ( 64), an expression for the critical amplitude in vertical direction can be extracted with However, it is not a necessary condition that the vertical acceleration amplitude has to be greater than the gravity constant g, see Sect. 4. With a sufficiently high coefficient of restitution, the vertical acceleration amplitude may also be smaller than g and the part is still conveyed, see Sect.5.1.
In Fig. 16, the comparison between the appearance of the second state of motion in the simulation model for several maximum initial heights z 0 and the prediction with the formula in Eq. ( 64) is shown for different Fig. 14 Multiple solutions of conveying system with different coefficients of restitution Fig. 15 Behavior of phases θ k of bouncing ball with different coefficients of restitution coefficients of restitution.It can be seen that the initial position z 0 ∈ [z 0,min , z 0,max ] is responsible for the appearance of the second state of motion.This initial position must be large enough so that the second state of motion appears at the predicted critical acceleration amplitude â * , see Fig. 16.In this case, the velocity of the point mass is equal to the velocity of the conveyor in feeding direction, see Fig. 14.If the initial position is too small, the real critical acceleration amplitude where the second state of motion appears is not equal to the predicted critical acceleration amplitude, see also Fig. 16.This means that the velocity at the first appearance of the second state of motion does not reach the velocity of the conveyor, see, e.g., Fig. 6.The phenomenon concerning the initial conditions and the corresponding state of motion will be studied later.In Fig. 14, it can be seen that the feeding velocity of the second state of motion is equal to the velocity of the Fig. 17 Conveying velocity of second state of motion depending on coefficient of restitution c r for several angles of throw ϕ conveyor in x-direction at the predicted critical acceleration amplitude â.Therefore, the feeding velocity of the second state of motion v * T can be expressed as With the predicted critical acceleration amplitude from Eq. ( 64), the feeding velocity of the part v T at the second state of motion can be calculated with which corresponds with the general approach in Eq. ( 28).
It is shown again that the feeding velocity v * T does not depend on the excitation amplitude.
In Fig. 17, the conveying velocity at the predicted critical acceleration amplitude in dependency of the coefficient of restitution c r and angle of throw ϕ is shown.The factor n is assumed to be 1, which means that the conveyor moves one period during the flight phase of the part.However, if n would be greater than 1, the feeding velocity increases by n times, see Eq. (67).In general, it can be seen that a lower coefficient of restitution results in a higher feeding velocity.Furthermore, the maximum possible feeding velocity increases if the angle of throw decreases.It has to be mentioned that the initial conditions have to be considered when applying this criterion.This is part of the next subsection.

Sensitivity of initial conditions
In this subsection, the sensitivity of the initial conditions of the conveying system is discussed.From the analyses in Fig. 16, it can be deduced that the real critical acceleration amplitude is also a function of the initial position.The predicted critical acceleration amplitude in Eq. (64) only describes the value above which the second state of motion may appear.However, the real critical acceleration amplitude ã * can be computed in general with where â * denotes the predicted critical acceleration amplitude from Eq. ( 64) and Δ â * denotes the unknown shift of the critical acceleration amplitude due to the initial conditions.
In [18], it is shown that the initial position z 0 is very sensitive and especially not independent because it has an influence on the time and the velocity at the first impact which are responsible for the steady state after the transient process.Therefore, the independent variables v 1 and ψ 1 are used in the following as initial conditions, which denote the velocity and the phase at the first impact.
In Fig. 18, two cases with different initial conditions v 1 and ψ 1 are shown.It is obvious that the real critical acceleration amplitude ã * depends on these, see Eq. ( 68).It can also be seen that the lower velocity (Sect.4) is present up to a higher excitation amplitude.Due to the initial conditions, sequences of microthrows and coupled normal motions can still occur even though the predicted critical amplitude has already been reached.
For a deeper understanding, especially of the shift of the critical acceleration amplitude, a variation of the independent initial conditions is performed.The real critical acceleration amplitude, where the second state of motion appears for the first time, is used as output variable.In Fig. 19, the results of this investigation are shown.
The white area in Fig. 19 shows that there is a large number of initial conditions where no critical acceleration amplitude exists, which means that the second state of motion does not appear.The colors in Fig. 19 indicate the critical acceleration amplitudes from â = 50 m/s 2 to â = 76 m/s 2 in several intervals.The contours between the transitions of the intervals show the same characteristic as the contours in [18] where the influence of the initial conditions is investigated at a constant acceleration amplitude â.The added value of the illustration in Fig. 19 is that the critical acceleration amplitude above which the second state Furthermore, the limits of the critical acceleration amplitude in Fig. 19 may also be estimated with the stability criterion from Eq. (40).In Fig. 20, the lower limit and upper limit for the respective coefficient of restitution are shown.The curve of the lower limit corresponds with the predicted critical acceleration amplitude in Fig. 16.The upper limit contains the new information about the maximum critical acceleration amplitude where the second state of motion appears for the first time.Between these two curves, the system shows a very sensitive but deterministic behavior, see Fig. 19.However, the sensitivity with respect to the initial conditions makes it very difficult to adjust the critical acceleration amplitude in practice.

Conclusion and outlook
Some approaches for the simulation and optimization of vibratory conveyor systems have already been presented in the past [9,11].In this work, a new approach for modeling and simulating vibratory conveyors is developed.The simulation model of the conveyor is based on the instantaneous bouncing ball model and allows to understand effects which appear in a conveying process.An advantage of the instantaneous formulation is that it leads to a massive reduction in computing time compared to a continuous formulation.
A special effect of vibratory conveyors is the occurrence of multiple states of motion at the same excitation amplitude [18].An analytical formula is developed to predict the first possible appearance of the second state of motion.This formula is verified with numerical experiments and the applicability in practice is shown.Furthermore, a formula is derived which allows to predict the conveying velocity in the second state of motion.However, the real critical excitation amplitude above which the second state of motion appears depends on the initial conditions.This is analyzed by numerical investigations in combination with a stability analysis of the corresponding bouncing ball model.
A further special effect is the microthrows which occur in a sequence of microthrows in general [4].The end of the sequence is reached if the relative velocity is equal to zero.Then the microthrow turns into a dragging process, for which the bouncing ball model has been extended.This dragging process is called coupled normal motion and is illustrated by an example, also with chaotic behavior.Again, analytical formulas are developed which allow an estimation of the duration of microthrows or coupled normal motions and the associated conveying velocities.These are also verified on the simulation example.
In future work, the developed formulas and statements will be applied to real measurements.These represent an enormous added value for practical applications.Furthermore, the impact model will be advanced and implemented in the multibody simulation tool HOTINT [31,32] which will result in a massive reduction of the computation time compared to the present continuous formulation.

Fig. 1 Fig. 2
Fig.1 Contact situation with point mass m on conveyor before and after the impact

Fig. 3
Fig. 3 Cases (a) and (b) of contact situation if v T < w T

Fig. 4
Fig. 4 Cases (c) and (d) of contact situation if v T > w T

Fig. 7
Fig. 7 Time signals of mean feeding velocity and z-position of part/conveyor of second state of motion

Fig. 8 Fig. 9
Fig. 8 Time signals of mean feeding velocity and z-position of part/conveyor in the first state of motion

Fig. 11 Fig. 12
Fig. 11 Feeding velocity and z-position of part/conveyor during coupled normal motion

Fig. 13
Fig. 13 Time signals of mean feeding velocity and z-position of part/conveyor in chaotic state of motion

Fig. 16
Fig.16 Comparison of analytic formula and simulation results with different initial positions z 0

Fig. 20
Fig. 20 Stability limits of critical acceleration amplitude of conveying process