Nonlinear and chaotic dynamics of a vibratory conveying system

In this work, a simulation model of a vibratory conveying system is presented. The simulation model is based on a continuous contact formulation in vertical direction which is extended by a friction force in horizontal direction to simulate a conveying process. In contrast with complex 3D simulation tools, it enables the understanding of previously unexplained phenomena such as multiple feeding velocities at the same excitation amplitude, which are observed in practical measurements. The parameters that have an influence on this effect are investigated, and a method for predicting and adjusting the occurrence of multiple solutions is developed. It is shown that the calibration of the system is very difficult in practice, as it depends significantly on the initial conditions which are difficult to reproduce and predict. It is also shown that the system can exhibit chaotic behavior in some configurations. These chaotic states are shown with the simulation model by means of parameter studies, and the point at which the system becomes chaotic is predicted with the method of Lyapunov exponents and fractal dimensions. Knowledge of the chaotic states can be used to calibrate the conveyor, as they depend only on the excitation and not on the initial conditions. The interdependencies of the initial conditions are also discussed in more detail. Therefore, this work provides a deeper understanding of complex conveying processes using a simple simulation model.


Introduction
In modern automation industry, vibratory conveying systems play a crucial role. Therefore, the performance of these systems becomes more and more important. As a main tool to improve these processes, dynamic simulation can be used. However, there are many open questions in conveying processes that are difficult to answer even with 3D simulation models. In particular, an effect has been observed in measurements where the parts are transported with different conveying speeds although the excitation of the conveyor is the same. Furthermore, trajectories are observed that suggest chaotic behavior. The major challenges addressed in this paper are the description of these effects and the interpretation of the results; especially, the transferability of the results into practical application is of high importance. Therefore, remarks about this transferability are given throughout different sections.
The design of linear feeding systems in industrial automation is mostly based on experimental results. However, this leads to a number of problems because many effects of the complex feeding process are not understood. Therefore, simplified calculation approaches were developed, e.g., in [1,2]. In [3], a numerical simulation model is presented which describes the feeding process with a discrete element method. A main challenge in modeling a conveying process is the contact between part and conveyor. In [4], a contact model for rigid bodies is developed which allows to model this contact with different types of damping and stiffness. The implementation of the contact model is based on the fundamentals of contact mechanics in [5,6]. The work of Czubak in [7] explores alternatives for a vibratory conveyor by utilizing a cost-effective electro-vibrator. This conveyor is able to efficiently transport small, loose materials at varying speeds and can also temporarily halt material flow without shutting down the drive. The analysis and simulation results of the conveyor's performance are encouraging and pave the way for additional laboratory testing. In [8], a new vibratory conveyor is investigated that focuses on precise material dosage. Its transport abilities in the resonance zone are analyzed using analytical and simulation methods. Furthermore, an optimal operating point is determined through simulation and verified with an industrial conveyor. Regarding the vertical component of a conveying process with a simple point mass, it is obvious that this has the same characteristic than the so-called bouncing ball problem. This problem is treated in [9][10][11][12], and many dynamic effects as period doubling or bifurcation are analyzed. More general investigations of nonlinear dynamic systems are presented in [13,14] and [15]. In [16], the nonlinear dynamic of ultrasonic transducers with impact contact is analyzed which is very close to the problem in this work.
A very complex field of nonlinear systems is the chaos research [17,18]. Chaotic states can also occur in mechanical systems, for example, in connection with liquid pumps [19,20]. In [21], the chaotic dynamics of repeated impacts in vibratory bowl feeders are treated. A stability analysis of a vibratory feeder and different states of motion is performed in [22]. A main challenge is to predict the states in which the system shows chaotic behavior. Therefore, many methods are developed, e.g., the Lyapunov exponent [23]. A further method to evaluate chaotic systems is fractal dimensions. These are known from the analysis of geometries, e.g., the Sierpinski triangle [24]. An extension to dynamical systems is presented in [25] where the change of the fractal dimension in time is evaluated. An overview of different methods to compute and evaluate fractal dimensions is found in [26,27]. In [28], methods for evaluating stability and chaos are analyzed in general. However, the importance of neural networks and machine learning has increased recently. In [29], an artificial neural network is trained that could distinguish chaotic and regular dynamics. The superiority over traditional methods and the robustness to varying control parameters is shown. In [30], the algebraic and spectral properties of horizontal visibility graphs associated with periodic and chaotic time series are investigated. In this work, measures for the quantification of the chaoticity of the process are investigated.
In Sect. 2, a reduced simulation model is presented that describes a point mass which is conveyed on a rigid plate. The contact between point mass and conveyor is modeled with a continuous contact formulation. The coupling with the horizontal direction for simulating the conveying process is represented with a friction force. A special feature of the simulation model is that the point mass switches between flight and contact phase. The solution of the non-smooth differential equation is performed using time integration.
A main effect observed in measurements of conveying processes is the occurrence of multiple feeding velocities at the same excitation amplitude. This effect is studied in Sect. 3. First, the occurrence of different states of motion is shown on a numerical example. Next, the simultaneous appearance of multiple solutions at the same excitation is reproduced. This effect is analyzed deeply, and a method is developed how to control the appearance of the respective state of motion by calibrating the initial conditions. This is a major challenge for practical applications because it is usually not possible to adjust the initial conditions as accurate as necessary.
A further challenge of studying the behavior of vibratory conveyors is the apparently chaotic behavior of the parts. In Sect. 4, methods are presented how to assess and proof the presence of chaotic states. The first method is the Lyapunov exponent which is discussed for discrete and continuous systems [25]. The second method is based on fractal dimensions which can be used to characterize special geometric shapes as the Cantor Set [31] or Sierpinsky triangle [32]. These methods are applied to the presented simulation model of the conveying system, and the chaotic behavior is assessed. Fig. 1 Transport of a 3D object in feeding direction (xcoordinate) by a vibratory conveyor moving in x-and z-direction However, the chaotic behavior has been observed in practice and 3D simulation already [33], but a proof is pending, which is part of Sect. 4. In addition, parameter sensitivity is investigated for assessing the transferability in practice.

A 2D-model for the simulation of linear feeding systems
In vibratory conveying systems, parts are transported on a vibrating channel as can be seen in Fig. 1. The conveyor N is excited in x-and z-direction which results in a movement of the part O in x-direction. In general, this conveying procedure can be applied to parts with complex 3D geometries. Due to increasing requirements, the simulation of such conveying systems becomes more and more important. In [4], a 3D simulation model was developed which allows the prediction of the behavior of the conveying parts. However, there appear some nonlinear effects which can be simulated but cannot be understood well using the complex simulation model. Therefore, a simplified model of a conveying process with a point mass is developed, see Fig. 2. This model should be able to describe the characteristics of effects appearing in measurements and 3D simulation.

Modeling of contact process in normal direction
In the following, the index p refers to the feeding part and the index s refers to the conveyor. The system distinguishes between two states, either the part is in contact (z p ≤ z s ) or not (z p > z s ). The mathematical distinction is modeled with the Heaviside-function θ (z p −z s ) which switches between 0 and 1 depending on the sign of the argument z p − z s . Assuming a gravity force with F g = −mg in z-direction, a stiffness force F k and a damping force F d , the differential equation for the conveying process can be written in general as The conveyor is generally excited with an acceleration amplitudeâ, with the angle of throw ϕ and a time-dependent cosine oscillation with the frequency f and a initial phase shift ψ 0 . The amplitude can be separated into the horizontal and vertical component withâ x =â cos(ϕ) andâ z =â sin(ϕ). Therefore, the full excitation in vertical direction reads The knowledge of the acceleration signalz s enables expressing the velocity and position signal witḣ and Next, the contact stiffness force F k and the contact damping force F d are investigated more closely. These forces are based on the implementation of the contact in the multibody simulation tool HOTINT, see [4].
The contact stiffness force F k is generally expressed with where k c denotes the stiffness. The exponent p enables to switch between Hooke's law ( p = 1) and Hertz' law ( p = 1.5), see [34].
The contact damping force F d is assumed as a velocity dependent force with where d c denotes the damping coefficient. This formulation is equivalent to a viscous damping, see [35]. The respective differential equation for the feeding part in vertical direction reads

Motion in translational direction
The modeling of the contact in Sect. 2.1 only considers the vertical movement of the part. However, the aim of a conveyor system is to move a part horizontally. Therefore, the vertical movement of the part from Sect. 2.1 has to be coupled with the horizontal direction. This coupling is realized with a friction coefficient μ which relates the contact force F c to the friction force F r , where The contact force F c in Eq. (8) is defined as the sum of stiffness and damping forces in the normal direction with It is not possible to evaluate Eq. (8) analytically because the differential equation Eq. (7) cannot be solved explicitly. As mentioned above, these differential equations are solved with a numerical time integration scheme. Therefore, the contact force F c is evaluated at every time step t i which enables to evaluate the norm of the friction force |F r | at every time step t i . For the consideration of the direction of the friction force, the velocity difference between part and conveyor in horizontal direction is used. The horizontal movement of the conveyor can be calculated analogously to Eq. (2) and (3) witḣ The sign of the friction force is equivalent to the sign of the velocity difference of conveyor and part in feeding direction. This results in Dividing the resulting friction force in Eq. (11) by the mass of the part m p results in the acceleration in feeding direction The numerical integration during the time integration process yields the feeding velocityẋ p which is a major quantity of a conveyor, see [36,37]. A further integration yields the position x p which is also suitable for some evaluations, especially sorting efficiency, see [38].

Analysis of multiple feeding velocities
In the following, investigations are performed with the presented simulation model. For the analyses, Hertz' law with p = 1.5 is used for modeling the stiffness force F k . The default set of excitation parameters is defined with the acceleration amplitudeâ = 40 m/s 2 , the angle of throw ϕ = 12 • , the frequency f = 50 Hz and the initial phase shift ψ 0 = 0 rad. The default set of contact parameters is defined with the stiffness coefficient k c = 5 · 10 5 N/m 1.5 , the damping coefficient d c = 0.3 Ns/m and the friction coefficient μ = 0.15. The mass of the part is m p = 0.7 · 10 −3 kg, and it starts at an initial position of z 0 = 0 mm with an initial velocityż 0 = 0 mm/s. The code for the simulation of the conveying process presented in Sect. 2 is implemented in the Python programming language. The equations are integrated with Euler's method and a step size of Δt = 5·10 −7 s which should guarantee convergent results. In the code, it is checked at each time step whether the part penetrates the conveyor. If this is the case, the contact algorithm is applied. If this is not the case, only gravity is applied and the part continues to fall down.

States of motion
In [33], it has been shown in measurement and simulation that difference states of motion of a part on a vibratory conveyor exist. The simplified model presented in Sect. 2 should represent this effect which enables a deeper consideration. First, parameter variations are performed to understand the principle behavior of the simulation model. As output variable, the mean feeding velocity is used. This mean value is calculated about In Fig. 3, the mean conveying velocity is computed for different acceleration amplitudesâ. It can be seen that the behavior changes significantly above a critical amplitudeâ = 50 m/s 2 where the second state of motion appears. This is equivalent to a vertical acceleration amplitude ofâ z = 50 sin(12 • ) = 10.40 m/s 2 which exceeds the gravitational acceleration.
In this work, the solutions below the critical acceleration amplitude are called first states of motion. The solutions above the critical acceleration amplitude are called second states of motion as long as the solutions are in a steady state and do not show chaotic behavior, which is analyzed in Sect. 4. Different acceleration amplitudesâ will also be called operating points.
In Fig. 4, the mean feeding velocity and the positions of conveyor and point mass in z-direction are shown, respectively, below and above the critical value. It is obvious that the feeding characteristic is changing between these two states of motion. Below the critical acceleration amplitude, the point mass moves periodically with the conveyor. Above the critical acceleration amplitude, the point mass has a jumping transport behavior.
In Fig. 5, the contact and friction forces of two states of motion are shown. It can be seen that the coupling between vertical and horizontal force results in changing the sign. This is caused by the relative velocities between part and conveyor, see Eq. (11).
For understanding different states of motion, an FFTtransformation of the velocity in vertical directionż p is applied. The results are shown in Fig. 6. It can be seen that multiple frequencies exist after the change to the higher velocity which is due to the hopping of the part.
In Fig. 7, the results of the same investigation with the velocity in horizontal direction are shown. It is shown that the first state of motion consists of more frequency components because the friction force is able to accelerate and decelerate the part. The second state of motion is dominated by a constant component because it is only a line with a few deflections during the contact process, see Fig. 4.
The differential equation in Eq. (7) leads to the assumption that the parameters mass m p , stiffness coefficient k c , damping coefficient d c and friction coefficient μ have an influence on the behavior of the system. Variations in the respective parameters have shown that they also influence the critical acceleration amplitude at which the system changes to another state of motion. However, these parameters are constant in practical application, which is why a variation in these is not analyzed further.
An interesting behavior is shown by a variation in the initial position z 0 . This parameter is very sensitive with respect to the output variable. A variation in the initial position at a constant acceleration amplitude may result in another state of motion. This means that several solutions exist at this operating point depending on the initial position. Furthermore, in practice, the initial position z 0 is the parameter that is most difficult to adjust accurately. Therefore, it is investigated more in detail in the following.

Results for different states of motion
In Sect. 3.1, it is shown that the simplified simulation model is able to represent different states of motion. A special effect discovered in [33] is that there occur multiple feeding velocities at the same operating point. It is assumed that these result from inaccuracies of the initial configurations which has been shown with the 3D simulation tool HOTINT in [33]. Therefore, the same investigations are performed with the simulation model, presented in Sect. 2. Simulations are performed at several operating points where the initial position is varied at every operating point in N steps. In Fig. 8, the results of the investigation are shown. The initial position was varied in N = 10 steps from z 0 = 0 mm to z 0 = 0.5 mm. It is obvious that in the range from 29 to 47 m/s 2 , multiple solutions exist. Fur-  thermore, the maximum velocity of the conveyor in x-direction is shown for comparing the feeding velocity of the part. However, it has to be mentioned that there appears chaotic behavior aboveâ ≈ 55m/s 2 . The chaotic behavior is studied in Sect. 4.
Next, the domain between 29 and 47 m/s 2 where two states of motion appear simultaneously is investigated in more detail. In Fig. 8, it can be seen that the variation in the initial position z 0 is responsible for the change of the state of motion at the respective acceleration amplitude. Figure 9 shows the dependency of the mean feeding velocity on the initial position z 0 for different acceleration amplitudesâ. It is obvious that the conveying speed increases abruptly at a critical ini-    Fig. 9 that the maximum critical initial position is 0.31 mm. If the initial position is larger than 0.31 mm, only the second state of motion occurs.
With regard to the occurring conveying speeds, it is also noticeable that the velocities of the first state of motion increase from 1 to 2 cm/s within an acceleration range of 18 m/s 2 . Looking at the velocities of the second state of motion, it can be seen that they are decreasing, although the accelerations are increasing. However, it must be noted that they only decrease by approximately 0.4 cm/s. It is also noticeable that the critical initial position does not depend linearly on the acceleration amplitudes, see Fig. 9. This is also obvious because of the horizontal distances between the curves which become smaller and smaller with increasing acceleration ampli-  tude. Therefore, the relation between the critical acceleration amplitudeâ and the respective critical initial position z 0 is investigated.
In Fig. 10, the critical initial position z 0 is plotted over the critical acceleration amplitudeâ. The nonlinear relationship which was determined with the interpretation of Fig. 9 is confirmed with this illustration.
The nonlinear relationship in Fig. 10 can be described with a function in the form where A, B and C are unknown parameters for the scaling of the regression curve. The parameters of the curve in Fig. 10 Fig. 10. The value of the error function is 4.14×10 −11 m 2 after the optimization. An essential question is how to choose the parameters such that the start and end of the respective state of motion in Fig. 8 can be controlled. In general, the initial condition z 0 is varied in an interval of the form From the previous investigation in Fig. 10, we know that a decreasing acceleration amplitude results in an increasing critical initial position. Therefore, it can be stated that an increase in the lower interval limit z L results in a lower critical acceleration amplitude which means that the first state of motion is no longer present at a lower acceleration amplitude.
With the upper interval limit z R , the acceleration amplitude, where the second state of motion appears, can be controlled. If z R increases, the second state of motion appears at a lower acceleration amplitude which means that the domain where two solutions exists has been extended.
For the verification of these assertions, a numerical example is performed where the end of the first and the start of the second state of motion should be modified. The requirement is that the first state of motion vanishes at 40 m/s 2 and the second state of motion begins at 35 m/s 2 . Therefore, the domain with multiple states of motion should be between 35 and 40 m/s 2 . To determine the interval limits of z 0 , Fig. 10 has to be evaluated. The end of the first state of motion should be at 40 m/s 2 . Therefore, the lower interval limit has to be at z 0 (40m/s 2 ) = 0.07 mm, see Fig. 10. The second state of motion should start at 35 m/s 2 , so the upper interval limit has to be at z 0 (35m/s 2 ) = 0.13 mm, see also Fig. 10.
The simulation with the initial conditions z 0 [0.07, 0.13] mm is shown in Fig. 11. It is obvious that the desired states of motion have occurred.
In this subsection, the nonlinear effect of multiple states of motion has been investigated and shown on different examples. The influence of the initial position z 0 is the focus of these analyses. This is due to the fact that z 0 has a practical relevance and has also been varied in [33]. However, there are also other initial conditions whose influences and interrelationships are investigated in the following.

Variation in initial conditions
In Fig. 12, the dependency of the states of motion on the initial position and vertical acceleration is shown. The colors distinguish between the first state of motion (blue), second state of motion (orange) and chaotic state of motion (green). Figure 12 summarizes the statements from Figs. 8, 9 and 10; especially, the frequency of the existence of solutions at lower and upper velocity between a = 28 m/s 2 and a = 47 m/s 2 (compare Fig. 9) and the nonlinear characteristic of it (compare Fig. 10) are illustrated. A novel information is that the initial position has no influence on the beginning of the chaotic states. Therefore, there exists a critical acceleration amplitude for switching to chaos which is independent of the initial position z 0 . For a generalization of the possible states of motion, the phase shift of the conveyor ψ 0 is varied additionally to the initial position z 0 . The initial velocityż 0 is set equal to zero. In Fig. 13, the results at an acceleration amplitude ofâ = 40 m/s 2 are shown.
The initial conditions of the analysis in Fig. 13 are the position z 0 , the velocityż 0 and the phase shift of the conveyor ψ 0 . It has to be mentioned that these initial conditions are not independent because a variation in the initial position z 0 has an influence on the time of the first impact t 1 and the velocity at the time of the first impact v 1 . Therefore, the initial position of the part is set equal to the initial position of the conveyor with Consequently, the initial velocityż 0 and the initial phase shift ψ 0 can be varied independently from each other and the whole phase space in 2D is covered. The results of the variation in the independent variablesż 0 . (16) With the point in time of the first impact t 1 , the corresponding velocitẏ z 1 = −gt 1 (17) and the corresponding phase of the conveyor can be computed. In order to use independent variables, z 1 and ψ 1 are set as initial conditionsż 0 and ψ 0 . Additionally, the initial position of the part z 0 has to be set equal to the initial position of the conveyor z s (0). In the investigations concerning the occurrence of multiple solutions, the initial position has been varied, see Fig. 8. However, the relations in Eqs. (17) and (18) make it possible to investigate which variations in the phase space result from the variations in the initial position z 0 . In Fig.15, it is shown that a variation in the initial position z 0 at a constant initial velocityż 0 and phase shift ψ 0 results in a straight line in the phase space. In the next section, the chaotic behavior of the conveying process is analyzed. With the criteria developed in Sect. 3.2, it is only possible to control the occurrence of the states of motion in a steady state, see Fig. 11. However, the investigations regarding the initial conditions in this subsection are significant for the analysis of chaotic behavior. This is the topic of the next section.

Analysis of chaotic behavior
The chaos theory refers to a subfield of nonlinear dynamics or dynamical systems that is not clearly defined. Essentially, it deals with special dynamical systems which evolution in time appears unpredictable, although the underlying equations are actually deterministic. This behavior is called deterministic chaos and arises when systems depend sensitively on the initial conditions [17,18]. Repetitions of an experiment can thereby lead to highly different measurement results due to minimal, hardly distinguishable initial conditions in the long-term behavior. However, chaos theory does not state that actually identical initial conditions would lead to different results.
An essential question is how to investigate and characterize dynamical chaotic systems. One possibility is to use the Poincaré map [39,40]. The Poincaré map is also used for describing the bouncing ball problem. Regarding the bouncing ball problem, a few typical chaotic phenomenons can be shown. One of it is the doubling of periods, investigated, e.g., in [41]. This effect describes that toward the chaotic range the oscillation period increases stepwise by a factor of two. Dif- ferent states in the bouncing ball problem are investigated in [9][10][11][12]. A special result of the investigations in [12] is that there exists a sufficient condition for the existence of a period one orbit.
A method to show the evidence of chaotic behavior of dynamical systems is the computation of Lyapunov exponents. One advantage of the method is that the indicator for the presence of chaotic behavior is only the sign of the Lyapunov exponent. A disadvantage is that the initial condition of the system has to be perturbed for getting a statement about the sensitivity [23]. A summary of the application of Lyapunov exponents can also be found in [42].
A further possibility to describe chaotic states is fractal dimensions. There exist several approaches to describe fractal dimensions. The basic idea is to reduce a set of n-dimensional points to a scalar value [43]. This is applied for example to images which show complex geometric behavior, e.g., the Cantor Set [31] or Sierpinski triangle [32]. In this work, it is used as an indicator for chaotic behavior.
In the following, the chaotic behavior of the conveyor is analyzed.

Observation of chaotic behavior of vibratory conveying systems
In Fig. 8, it can be seen that the simulation model shows a chaotic behavior if the acceleration amplitude is abovê a ≈ 60 m/s 2 . In Fig. 16, the time signals of a simulation with an acceleration amplitude ofâ = 80 m/s 2 are shown. Compared to the investigations in Fig. 4, it is obvious that the characteristic of the trajectory is equivalent to a state of the second state of motion. However, even with a longer simulation time, no steady state is achieved which is the first indication of chaos. Furthermore, in Fig. 16 the associated contact and friction force are shown. It can also be recognized that no periodic behavior is evident. Next, the sensitivity of the solution to the initial condition is investigated. In Fig. 17, the dependency of the mean velocity on the initial position z 0 is shown for several acceleration amplitudesâ i . It is obvious that chaotic behavior occurs, especially if it is compared with Fig. 9.
In Fig. 18, the possible states of motion are visualized over the phase shift ψ and the initial position z 0 . It can be seen that the contours are similar to Fig. 13 with the difference that chaotic states may appear. Due to the dependency of the initial position z 0 , the independent parameters ψ 0 andż 0 are chosen for the variation. The necessary condition that the initial position of the conveyor z 0 has to be equal to the initial position of the conveyor z s (0) is considered. Comparing Figs. 19 to 14, it is obvious that the characteristic is again very similar. The difference is the appearance of chaos in Fig 19. Furthermore, it can be seen in Fig. 19 that the phase shift may significantly determine which state the system reaches. The results from Fig. 18 could be converted with the relations described in Eq. (17) and (18) to the representation in Fig. 19. It is not possible to convert Fig. 19 into the representation in Fig. 18 because in Fig. 18 it is assumed that the initial velocityż 0 is equal to zero. Therefore, a conversion from x 0 = (ż 0 , ψ 0 ) to x 0 = (z 0 , ψ 0 ) is only meaningful for some configurations. However, the evidence of chaotic behavior is still pending. This is performed with Lyapunov exponents and fractal dimensions in the following.

Evidence of chaotic behavior
A method to assess the existence of chaotic behavior of dynamical systems is the method of Lyapunov exponents. The Lyapunov exponent for a discrete dynamical system x n+1 = f (x n ) is defined as and describes the exponential separation of two nearby trajectories originally separated by distance ε [42]. For a continuous system, the definition of the Lyapunov exponent in Eq. (19) has to be extended. In general, a continuous dynamical system is written aṡ with x(0) = x 0 andẋ(0) =ẋ 0 . Explicitly considering the initial condition x 0 and assumingẋ 0 = 0, the differential equation in Eq. (20) has the solution x(t, x 0 ). Therefore, the Lyapunov exponent from Eq. (19) is redefined as The derivative term in Eq. (21) may be written as Therefore, the Lyapunov exponent for continuous dynamical systems is defined as For a numerical calculation of the Lyapunov exponent of continuous dynamical systems, the perturbation in the initial conditions ε has to be chosen very small for evaluating Eq. (23). Furthermore, the run-up process of the transient procedure is skipped for the computation of the Lyapunov exponents. This is ensured by using only the last 25 % of the measuring distance, compare Sect. 3. First, the Lyapunov exponent of the system is calculated with Eq. (21), where x(t, x 0 ) is chosen as the position of the point mass in z-direction z p (t, z 0 ). For this analysis, the data from Fig. 8 are used to compute the numerical derivative ∂z p (t,z 0 ) ∂z 0 at z 0 = 0. With this approach, the perturbation is equal to where Δz 0 is the interval of initial conditions and N z is the number of steps in this interval. The resulting Lyapunov exponent, evaluated at several acceleration amplitudesâ i , is shown in Fig. 20. It is recognizable that the Lyapunov exponent becomes greater than zero atâ = 42 m/s 2 . We know from Fig. 8 that the system does not show a chaotic behavior at this operating point but only the occurrence of two states of motion. Betweenâ = 50 andâ = 60 m/s 2 , the Lyapunov exponent changes the sign repeatedly, although the behavior is already chaotic, see Fig. 8. Therefore, the assessment of chaos with the already existing data is not enough. Furthermore, it has to be mentioned that the initial position z 0 of the system in Fig. 20 is used as perturbation parameter. As shown in Fig. 15, the initial condition is not an independent variable wherefore it should not be used as perturbation parameter. Therefore, the independent phase shift of the conveyor ψ 0 is chosen as perturbation parameter. The respective initial velocity is chosen as a constant value and the initial position of the point mass is equal to the initial position of the conveyor.
Next, the Lyapunov exponent of the time continuous system is calculated with a small perturbation ε in the phase shift ψ 0 with the forward difference quotient in Eq. (23). In Fig. 21, the influence of the perturbation ε on the time signals of the measures in horizontal direction is shown.
The results of the calculation of the Lyapunov exponent of the output variableẋ p are shown in Fig. 22 for several acceleration amplitudesâ i . It is obvious that the Lyapunov exponent is greater zero above a = 50 m/s 2 which is the assessment criterion for chaotic behavior. This statement is consistent with the observations in Fig. 8. Therefore, chaos has been demonstrated using the presented method. The reason for the deviations in Fig. 20 is the size of the perturbation ε and the chosen perturbation variable z 0 which is not independent as shown in Sect. 3.3.
After verifying the chaotic behavior with the Lyapunov exponent, the method of fractal dimensions is presented. In general, the fractal dimension of a set is a generalization of the dimension concept of geometric objects such as curves and surfaces. The special property is that the fractal dimension does not have to be an integer. In this work, the fractal dimension is described as Minkowski-Bouligand dimension, also called boxcounting dimension [44].
For calculating the Minkowski-Bouligand dimension D 2 for a set E given by the points {X i |i = 1, ..., N }, a correlation integral C(R) has to be defined. This correlation integral describes the number of sets C which have a smaller distance to each other than R. The evaluation of C(R) is done with grids of space, which is where the name box-counting comes from. The derivative of the correlation integral C(R) corresponds to the fractal dimension D 2 [25] which is defined as In Fig. 23, the phase portraits at the three states of motion (lower velocity, upper velocity, chaos) are shown. It has to be mentioned that the vertical position z p (t) and vertical velocity v z (t) are used for the phase portrait. The reason is that the position in horizontal direction x(t) increases with time because of the feeding process (see Fig. 1). Therefore, a closed trajectory is not possible. The time t = 1 s has been chosen as the  Fig. 23. It can be seen in Fig. 23 that the phase portraits of the states with lower (left) and upper velocity (middle) have closed contours which means that they should have the dimension 1. The phase portrait of the chaotic states does not consist of closed contours. Therefore, the longer the simulation time, the more the space is filled and the fractal dimension should be between 1 and 2.
In Fig. 24, the fractal dimension of the phase portrait at several acceleration amplitudes is shown. It is obvious that the beginning of the chaotic behavior can be predicted with this method, compare Fig. 22. It can be seen that the fractal dimension fluctuates a lot. This is also a disadvantage in comparison to the Lyapunov exponent, since with the latter the sign is sufficient for an assessment. However, for assessing the beginning of the chaotic behavior, this method is suitable.

Conclusion and outlook
The simulation of vibratory conveying systems is a main challenge, especially due to the nonlinear effects. In [4] and [8], some approaches for simulating and optimizing vibratory conveying systems are presented. A special challenge here is the correct modeling of various effects. In this work, the dynamic behavior of a vibratory conveying system is analyzed in general. The presented simulation model of a vibratory feeder allows a deeper understanding of the effects which occur in a conveying process. The first effect is the occurrence of multiple feeding velocities at the same operating point. The dependencies of the model parameters are investigated and a procedure is presented which allows to calibrate the appearance of the respective state of motion. It is shown that the initial conditions are very sensitive and are mainly responsible for this effect. However, it is very difficult to use this in practice because the initial conditions cannot be reproduced and adjusted in the measurement as accurately as necessary. The second effect is the occurrence of chaotic behavior if the excitation exceeds a critical value. The chaotic behavior is observed in the simulation results and is proved with suitable methods, the Lyapunov exponent [16] and fractal dimensions [25]. A further result is that the critical point at which the system changes to chaos is mainly depending on the excitation amplitude. This enables the adjustment of the vibratory feeder in practice.
In future work, the understanding of the sensitivity of the states of motion is to be expanded. For this purpose, analytical approaches based on simple models will also be developed. The results from these approaches will be verified in practice together with the results from this work. It has to be mentioned that the presented simulation model contains only a point mass which is an obvious limitation of the presented work. However, if rigid bodies would be used, there may be additional effects that we can not investigate with this model. Therefore, the developed methods will be transferred to the multibody simulation tool HOTINT [45] and used for parameter identification of vibratory conveyors [33]. This is intended to generalize the knowledge gained.
Funding Open access funding provided by Johannes Kepler University Linz. This work has been supported by the COMET-K2 Center of the Linz Center of Mechatronics (LCM) funded by the Austrian federal government and the federal state of Upper Austria. The funder name is FFG (Österreichische Forschungsförderungsgesellschaft) and the Grant ID is 886468.

Data availability statement
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors have no competing interests to declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.