Development of a mathematical model of a vibrating polyfrequency screen as a dynamic system with distributed parameters

. A mathematical model of a vibrating polyfrequency screen as a dynamic system with distributed parameters has been developed. The dynamic system of solids of finite sizes was chosen as the design scheme for the screen: framework, sieves with bulk material and impactors, the contact interaction of which occurs through two-side bonds and collisions on surface areas that have elastic-damping coverings. It is shown that a change in the amplitude of the exciting force has a significant effect on the dynamics of vibration impactors of a polyfrequency vibrating. There is an amplitude value at which the impactor passes from the mode without interacting with elastic bonds to the vibro-impact mode. The impactor movements begin to change disproportionately altered by the exciting force amplitude. It is shown that the start of the impactors in the screen substantially depends on the exciting force. Changes in the amplitude of the exciting force make it possible to achieve chaotic oscillations of impactors, which in turn leads to oscillations of screen surfaces with a continuous frequency spectrum, i.e. to the operation mode of the screen, which is most appropriate for dehydration and separation of fine mineral fractions.


Introduction
Efficient extraction, dehydration and enrichment of fine fractions of mineral raw materials is possible on vibrating screens, the screening surface of which make polyfrequency oscillations with a continuous infinite frequency spectrum and with accelerations up to thousands of m/s 2 and the necessary level of acceleration of the screening surface have not been fully performed. Therefore, modeling and setting parameters of a vibrating polyfrequency screen for effective separation by size and dehydration of small fractions of mineral raw materials is relevant.
In [2], for the research of nonlinear polyfrequency vibrations of vibrating screens, a model consisting of concentrated bodies was developed. Each body could only perform vertical vibrations, i.e. it had one degree of freedom. Computer modeling of the dynamic behavior of a discrete system made it possible to establish a number of important laws of nonlinear polyfrequency vibrations of vibration screens with one impactor [3]. However, if there are two or more impactors in the screen, the model consisting of concentrated bodies does not reflect the essential structural design feature of the screen -the distributed nature of the contact interaction between the impactors and the sieve. Therefore, to research the dynamics of vibrating screens with several impactors, it is necessary to develop a mathematical model of the screen as a distributed parameters system (Fig. 1).
The most common mathematical model of a mechanical system with distributed parameters is a system of contacting three-dimensional deformable bodies with non-linear conditions of contact interaction, taking into account the one-sided character of the bonds, the presence of friction and coatings of contacting surfaces, as well as other tribological factors [4,5]. Various rheological models of deformable bodies can be used to describe the processes of energy dissipation during vibrations. The most common are the models of viscoelastic bodies by Voigt, Maxwell, and Kelvin (the standard model of a viscoelastic body) [6]. In [7], a generalized nonlinear model is proposed for taking into account energy dissipation during vibrations on the basis of an invariant representation of hysteresis loops, and in [8] an approach using a complex module is described.
In the general, the equations of mechanical systems dynamics with distributed parameters are systems of partial differential equations and have an infinite number of degrees of freedom [9]. It is very difficult to obtain an analytical solution of dynamic contact tasks for such systems of equations due to the significant nonlinearity of the boundary conditions; therefore, numerical methods are usually used to solve them [10][11][12].
At the first stage of formation a numerical solution, it is necessary to perform a halfsampling of the task by spatial coordinates, i.e. to approximate a distributed parameters system by a discrete model with a finite number of degrees of freedom. In this case, the original task is reduced to a system of nonlinear ordinary differential equations.
To solve the task obtained as a result of half-sampling by the spatial coordinates of the nonlinear ordinary differential equations system, as a rule, various numerical integration algorithms are used [13][14][15][16]. However, the large dimension of these systems of equations requires powerful computational resources to solve them, in particular, the use of parallel computing systems [17].
When conducting applied research, a rational simplification of mathematical models is possible, which does not lead to a loss of the results reliability of quantitative analysis, but makes these results more visible and interpretable in practical terms. Therefore, when solving many applied problems arising in engineering practice, it is advisable to use simplified mathematical models that take into account the shape, geometric dimensions, stiffness ratios, and ranges of other parameters variation of the such distributed parameters system objects.
One of the main methods of simplifying mathematical models is associated with the adoption and justification of the so-called kinematic hypotheses that establish certain restrictions on displacements, which, combined with the use of the principle of possible displacements and analysis of the possible work of the applied external load, makes it possible to identify in the calculated scheme of the considered technical object and its mathematical model the most significant factors and neglect the influence of less significant. The result of this simplification is the design schemes of the rod, membrane, plate, and shell that are widely used in applied mechanics [18].
η n ζ n -coordinate system are tightly bound with the bodies; m 1 , m 2 , m 3 , m n -bodies masses, kg; c, c p101 , c p103 , c p104 -stiffness coefficients of the two-side bonds, N/m; b, b p101 , b p103 , b p104 -damping coefficients of the two-side bonds, kg/s; δ 211 , δ 214 , δ n11 , δ 121 , δ 124 , δ 1n1 , δ 131 , δ 134 , δ 1n3 , δ 1n4 , δ n14 , δ n13 -initial gaps, m; l x , l y -length of one-side bound along x and y axis, respectively, m; Q -amplitude of harmonic force, m; ω -frequency of harmonic force, rad/s; t -time, s; y q Another common technique for simplifying mathematical models is used for structures whose elements have significantly different stiffnesses. In this case, one part of the structural elements is modeled by absolutely solid (stiff) bodies, and the other part is modeled by deformable bodies. Such an approach is used, for example, in researches of the dynamics of vibrating machines with elastic elements [19].
-displacement in the application of force relative to the origin.

Method
In developing the mathematical model of the screen, the following assumptions were made: -taking into account that the structural elements of the screen have substantially different stiffness, one part of the structural elements is modeled by absolutely solid bodies, and the other part is by Voigt's elastic-damping bodies; -a framework, impactors, as well as a sieve with bulk material are modeled by absolutely solid (stiff) bodies of finite sizes; -contact interaction of solids bodies occurs through holding (two-side) and non-holding (one-side) bonds; -two-side bonds are modeled by non-linear elastic-damping Voigt's elements, which are attached to a pair of solids at particular points on their surface. In the initial state of the system, the lines of action of all two-side bonds are directed vertically; -one-side bonds of solids are ideal non-holding bonds, i.e. in the case of collisions of solids, the reactions are directed along the normals to the contact surfaces, there are no friction forces; -the shape and size of the limiting contact zones on the surface of solids, through which one-way contact interaction (collision) with other solids occurs, can be indicated from geometric considerations. Actual areas of contact may vary during movement; -areas of the surface of solids along which one-way contact interaction (collision) occurs, in the initial state of the system are parallel to the horizontal plane, i.e. the normals to these surfaces are directed vertically; -the contacting sections of the surfaces of solids have surface coatings that are modeled by thin viscoelastic Winkler-type layers, i.e. only the deformation of the transverse compression of the coating along the normal direction to the surface of the solid is taken into account, and the viscoelastic properties of the coatings are described by the Voigt model; -the initial gaps between the contacting sections of the surfaces of solids, taking into account the coatings, are small in comparison with the linear dimensions of the bodies; -the value of the transverse compression of the surface coating along the normal direction to the surface of the solid is small compared with the linear dimensions of the body; -the vibrations of the screen are excited by inertial oscillators, generating vertically directed constrain forces that change in time according to a harmonic law. Oscillators have unlimited power and are attached to solids at particular points on their surface; -ideal bonds exclude horizontal movements of solids and their rotation relative to the vertical axis; -for each solid, one of the main axes of inertia is directed vertically. Taking into account the assumptions made, the dynamic system N≥3 of solids of finite sizes -a framework (body S 1 ), a sieve with bulk material (body S 2 ), and impactors (body ), whose contact interaction occurs through two-side bonds (Voigt nonlinear elastic-damping elements) and collisions on surface areas with elastic-damping coatings is selected as the calculation scheme for the screen.
The position of the bodies of the dynamical system will be determined in a fixed Cartesian coordinate system Oxyz, hereinafter referred to as the main system. We assume that the axis Oz is directed vertically upwards, and the axes Ox and Oy are located in the horizontal plane. We also introduce for each solid S n N n , 1 = , , a moving coordinate system O n ξ n η n ζ n having a origin (pole) O n in the center of mass of the body. The axes of the coordinate system O n ξ n η n ζ n are tightly bound with the body and are its main axes of inertia. Taking into account the assumptions made above, we assume that the coordinate axes O n ζ n of all bodies S n N n , 1 = , are directed vertically upward. All of the above coordinate systems are right-hand.
In the general case, the position of the moving coordinate system O n ξ n η n ζ n n on OO r =  relative to the main coordinate system Oxyz is determined by the vector and the three angles of rotation that are made in a certain sequence around the appropriately selected axes. The angles of rotation are called Euler angles [20]. For a given sequence of rotations, the angles are determined definitely. In the general case, finite rotations are noncommutative. In the particular case of small transformations, rotations have the commutativity property, i.e. for a small transformation, the values of small rotation angles do not depend on their sequence [21].
Further, let us agree with the positive direction of rotation to call the rotation carried out in the right-hand coordinate system counterclockwise for the observer looking from the positive tip of the axis around which the rotation is made.
In accordance with the assumption made above about the existence of ideal two-side bonds, excluding horizontal movements of solids and their rotation about the vertical axis, we assume that the centers of mass (poles) O n of all bodies S n N n , 1 = , , of the dynamic system can make small movements only along the axis Oz of the main system Oxyz, there are no displacements along the axes Ox and Oy, and the bodies S n themselves can make small rotations with relative to the axes O n ξ n and O n η n , at the same time, rotation of the bodies relative to the axes O n ζ n is excluded. Therefore, each solid body S n N n , 1 = , , of the considered dynamic system has three degrees of freedom: n w -movement of the center of mass (pole) O n along the axis Oz; α n -angle of rotation relative to the axis O n ξ n ; β n -angle of rotation relative to the axis O n η n Taking into account the accepted assumptions regarding the kinematics of the motion of bodies, the components of the displacement vector of an casual point of a solid body S where u, v, w -displacement vector components in the main system Oxyz, m; ξ n , η n , ζ ncoordinates of a point in a coordinate system O n ξ n η n ζ n tightly bound with the body, m; ϑ n -angle of rotation of the coordinate system O n ξ n η n ζ n Consequently, for solids of finite dimensions, the displacements of all their points will be small. relative to the coordinate axis Oz.
Thus, the position of the considered dynamic system is determined by the 3N generalized coordinates (degrees of freedom): Let us agree on the following designation for generalized coordinates: , i.e. the displacements of the centers of mass of solids along the vertical axis and the angles of their rotations relative to the horizontal axes are selected as generalized coordinates.
and form a coordinate vector We obtain the equations of motion of the dynamical system in the form of second-order Lagrange equations, which have the form [25] where E k -kinetic energy of the system, J; E p -system potential energy, J; q i i q  -generalized coordinates, m; -generalized speeds, m/s; Ф -dissipative function; F i -generalized non-potential force related to the generalized coordinate q i We present the reaction P , N.
where c pn , d pn , λ pn -stiffness characteristics coefficients of the n-th two-side bond, N/m; b pn -damping coefficient of the n-th two-side bond, kg/s; e pn pn e  -deformation of the n-th two-side bond, m; -deformation rate of the n-th two-side bond, m/s.

Moments of reaction P pn
relative to the coordinate axes are calculated by the formulas: where M pnξ1 -reaction moment P pn relative to the axis O i ξ i , i=I p1 (n), N⋅m; M pnη1reaction moment P pn relative to the axis O i η i , i=I p1 (n), N⋅m; M pnξ2 -reaction moment P pn relative to the axis O j ξ j , j=I p2 (n), N⋅m; M pnη2 -reaction moment P pn relative to the axis , j=I p2 n (n) -numbers of a pair of solids contact with -th two-side bond; I p1 (n), I p2 In accordance with the above assumptions, the shape and size of the limiting contact zones on the solids surface, according to which there is one-side contact interaction (collision) with other solids, are indicated from geometric considerations. We emphasize that only the limiting sizes of the contact zones are assumed to be known, the actual contact areas change during movement and must be determined when modeling the dynamic behavior of the system.
(n) -the integral functions of the integer argument, whose values are equal to the numbers of the pair of solids on which the n-th two-side bond is superimposed.
We also introduce two integral functions of the integer argument I r1 (n) and I r2 , whose values are equal to the numbers of a pair of solids contact with the n-th one-side bond.
The part of the body surface S i , i=I r1 (n), on which its contact with the body S j , j=I r2 (n) is possible, is denoted Ω n1 . Similarly, the part of the surface of the body S j along which its contact with the body S i is possible, is denoted Ω n2 . It is assumed that the external normal to the surface Ω n1 in the initial equilibrium state of the system is directed vertically downward, and the external normal to the surface Ω n2 It was shown in [22] that, at small displacements, to a first approximation, the form of the deformed body boundary is determined by the normal displacements of the points lying on it. Let the points Z -vertically upward. rn1 ∈Ω n1 and Z rn2 ∈Ω n2 have the same coordinates x and y in the main coordinate system Oxyz in the initial state of the system. Then the sum of the transverse compressions of surface coatings at these points of contacting surfaces Ω n1 and Ω n2 the n-th one-side connection can be calculated by the formula , where s rn -sum of transverse compressions of surface coatings, m; δ rn -initial gap between points Z rn1 and Z rn2 along the normal to the surface Ω n1 , m; w(Z rn1 ), w(Z rn2 )vertical movement of points Z rn1 and Z rn2 respectively, m; H(⋅) -Heaviside function (13) Note that the presence of the Heaviside function in expression (12) ensures that the condition of non-positivity of the value of the transverse compression of the coatings of the touching surfaces of a one-sided (non-holding) bond is satisfied: The rate of change of the transverse compression of the surface coatings of a one-side bond is calculated by the formula where rn s  -rate of change of transverse compression of surface coatings of n-th one-side bond, m/s; -vertical movement of points Z rn1 and Z rn2 Substituting expressions for the vertical displacement components (3) in (12) and (15), we obtain respectively, m. where where σ rn -reactive pressure, N/m 2 Moments of reaction P where M rnξ1 -reaction moment P rn relative to the axis O i ξ i , i=I r1 (n), N⋅m; M rnη1 -reaction moment P rn relative to the axis O i η i , i=I r1 (n), N⋅m; M rnξ2 -reaction moment P rn relative to the axis O j ξ j , j=I r2 (n), N⋅m; M rnη2 -reaction moment P rn relative to the axis O j η j , j=I r2 The sum of the elementary works of forces applied to a solid S (n), N⋅m.
where δΨ n -sum of elementary works of forces, J; N fn -number of inertial oscillators attached to a solid S n ; f ni -amplitude of i-th harmonic force, m; ω ni -frequency of i-th harmonic force, rad/s; ϕ ni -initial phase of i -th harmonic force, rad; t -time, s; δw(C ni )virtual vertical movement of a point C ni , in which the i-th harmonic force is enforced, m.

Substituting in (24) the expression for the vertical component of displacements (3), we obtain
where ξ fni and η fni -point Z fni From (25), taking into account (4), it follows that coordinates, in which the i-th harmonic force is enforced, m.
where F 3n-2 , F 3n-1 , F 3n -generalized non-potential forces referred to generalized coordinates q 3n-2 , q 3n-2 , q 3n N n , 1 = , respectively, N; . As a result, the system of equations of motion of the dynamic system takes the form: where m i -body mass, kg; J iξ -moment of inertia of the body relative to the axis O i ξ, kg/m 2 ; J iη -moment of inertia of the body relative to the axis O i η, kg/m 2 ; P pn -reactions of two-side bonds calculated by the formula (7), N; P rn -reactions of one-side bonds calculated by the formula (19), N; L p1 (i), L p2 (i) -sets of numbers of two-side bonds superimposed on a solid S L r1 (i), L r2 (i) -sets of numbers of one-side bonds into which a solid enters S If in (25)-(26) any of the sets L p1 (i), L p2 (i), L r1 (i), L r2 , is empty, then the appropriate sum is assumed to be zero. The system of differential equations (29)-(31) must be complete with initial conditions. If at the initial moment of time t=0 the dynamic system was at rest, then the initial conditions have the form: where . It is very difficult to obtain an analytical solution to the Cauchy problem (29)-(31), (36)-(38) due to its significant nonlinearity, therefore, it is advisable to use numerical integration algorithms to solve it. The computational time integration algorithm for the equations of motion (29)-(31) is based on the use of three-layer difference schemes with weights [21].
The use of a computational algorithm for integrating the equations of motion of the dynamical system allows one to obtain time series {Q k } describing the generalized coordinates of the dynamical system at discrete time moments t k The developed mathematical model of a screen with distributed parameters and a computational algorithm for its numerical research are implemented in the form of a computer program.
, which, as a rule, are taken at equal time intervals h, called the sampling period.

Research results and discussion
The researches of a model of a polyfrequency vibrating screen with a single impactor were carried out.
In Figure 2 shows a diagram of a vibrating polyfrequency screen with one impactor. The screen consists of a framework with mass m 1 , which is attached to a fixed base by means of support bonds with stiffness c p10 and viscosity b p10 . A sieve of mass m 2 with stiffness c p12 and viscosity b p12 is fixed on the framework. An impactor with mass m 3 is also attached to the framework using elastic bonds with stiffness c p31 and viscosity b p31 . Also, elastic limiters with stiffnesses c r31 , c r13 and viscosities b r31 , b r13 , respectively, are installed above and below the impactor. The gap between the upper limiter and the impactor -δ 13 , between the lower limiter and the impactor -δ 31 . Stiffness of the lining of the impactor -c r23 , viscosity of the lining -b r23 . The gap between the sieve and the impactor -δ 23 . The screen is set in motion by an external force with amplitude q 1 and frequency of excitation ω 1 In the Table 1 shows the parameters of the screen at which researches were performed.
. Table 1 1000 , kg/s The gap between the impactor and the sieve δ 23 0 , m Excitation force frequency ω 1 156 , rad/s In Figure 3 shows the dependence of the maximum displacements of the impactor on the amplitude of the exciting force q 1 Initially, with an increase in the amplitude of the exciting force from 0 to 48000 N, the amplitude of the displacements also gradually increases, while the vibration mode is unshocked. However, at a value of 49000 N, there is a sharp increase in the movements of the impactor and the transition of the impactor to the vibro-impact mode. This mode is characterized by a sharp change in the amplitudes of the impactor's displacements depending on the exciting force at q . Gray color on the graph shows the upper elastic limiter. Thus, when the curve of maximum displacements of the mass is under a straight line corresponding to the position of the limiter, i.e. the movement of the impactor is less than the gap, the vibration mode is unshocked, but if the curve of maximum displacements above the limiter, this indicates a vibro-shock mode of the impactor operation. ≥ 49000 N. In contrast to a single-mass system with elastic limiters, in the case of a single-impactor screen, the value of the force amplitude increases at which the vibro-impact mode begins, and there is also no area in which the regime is vibro-impact, but the displacements increase in proportion to the increase in the exciting force. In other words, =49000 N. Moreover , the difference between II q 1 ′ and I q 1 ′ ′ is relatively small (49000 and 52000 N, respectively). The researches of a model of a polyfrequency vibrating screen with two impactors were made. The above results of the researches made it possible to establish that the influence of the amplitude of the exciting force on the vibration mode of the impactor, and, consequently, the whole screen, is significant. The next step is the researches of a model of a vibrating polyfrequency screen with two impactors to establish whether there is a difference in the dynamics of this system and the one described above.
at which the impactor from the mode without interacting with the elastic limiters passes to the vibro-impact mode. In this case, the movements of the impactor begin to change disproportionately to the change in the amplitude of the exciting force.
The scheme of a vibrating polyfrequency screen with two impactors is similar to the scheme shown in Figure 2. The screen consists of a framework with mass m 1 , which is attached to a fixed base by means of support bonds with stiffness c p10 and viscosity b p10 . A sieve with the mass m 2 , stiffness c p12 and viscosity b p12 is fixed on the framework. Impactors with masses m 3 and m 4 are also attached to the framework using elastic bonds with stiffnesses c p31 , c p41 and viscosities b p31 , b p41 . Also, elastic limiters with stiffnesses c r31 , c r13 , c r41 , c r14 and viscosities b r31 , b r13 , b r41 , b r14 , respectively, are installed above and below the impactors. The gaps between the upper limiters and the impactors -δ 13 , δ 14 , between the lower limiters and the impactors -δ 31 , δ 41 . The stiffness of the lining of the impactors -c r23 and c r24 , the viscosities of the lining -b r23 and b r24 . Gaps between sieve and impactors δ 23 and δ 24 . The screen is set in motion by an external force with amplitude q 1 and frequency of excitation ω 1 . In the Table 2 shows the parameters of the screen at which researches were performed. To facilitate the task, the parameters of impactors were taken equal.
In Figure 4 shows the dependence of the maximum displacements of the impactors on the amplitude of the exciting force q 1  The upper elastic limiters are displayed in gray on the graph. Thus, when the curves of the maximum displacements of the impactors are under a straight line corresponding to the position of the limiter, i.e. the displacements of the impactors are less than the gap, the vibration mode is unshocked, but if the curves of the maximum displacements above the limiter, this indicates a vibro-impact mode of impactors operation.
From the graph shown in Figure 4, it can be seen that for a small value of the exciting force, the difference in the vibrations of both impactors is insignificant. However, with force increasing, this difference also increases. In this case, up to the value q 1 =37000 N, both impactors oscillate in unshocked mode. With force values greater than 37000 N, one of the impactors enters into the vibro-impact mode, and the second stay in unshocked. It should be noted that the impactor, operating in vibro-impact mode, is constantly changing. With a further increase in the amplitude of the exciting force, modes appear in which both impactors interact with elastic limiters, and the more q 1 In Figures 5-7 show phase diagrams of vibrations of impactors for various values of the amplitude of the exciting force.
, the more often these modes occur.
From the diagrams given in Figure 5 it can be seen that both impactors oscillate exactly the same in unshocked mode with a period equal to 2T (where T -time of oscillation period, s). The oscillation mode to which the phase diagrams in Figure 6 is non-periodic for both impactors. The Poincare sections shown in Figure 8, indicate that this is a regime of deterministic chaos. Such modes are characteristic for a wide range of values q =5000 N. 1 In Figure 7 shows the phase diagrams of the 4T-periodic regime, which is characteristic of both impactors. This is a vibro-impact mode, characteristic for values of q . 1 Thus, we can conclude that the start of both impactors in a vibrating polyfrequency screen significantly depends on the exciting force, i.e. the problem with the absence of a vibro-impact mode in one or both of the impactors can be solved by increasing the amplitude of the exciting force. In addition, by changing the amplitude of the exciting force, it is possible to achieve chaotic modes of impactors vibration, which in turn leads to vibrations of sieve surfaces with a continuous frequency spectrum, that is, to modes of screen operation, most conductive for dehydration and separation of fine fractions of minerals [1][2][3][23][24][25].   1. A mathematical model of a vibrating polyfrequency screen as a dynamic system with distributed parameters has been developed, which takes into account the presence of two or more impactors in the screen and reflects the significant structural feature of the screen -the distributed nature of the contact interaction between the impactors and the sieve. A dynamic system of finite-sized solids -a framework, a sieve with bulk material, and impactors, the contact interaction of which occurs through two-side bonds (non-linear elastic-damping Voigt elements) and collisions on surface areas with elastic-damping coatings -was chosen as the calculation scheme for the screen. The use of a computational algorithm for integrating the equations of motion of the dynamical system allows to obtain time series describing the generalized coordinates of the dynamical system at discrete time moments, which, as a rule, are taken at equal time intervals, called the sampling period. The developed mathematical model of a screen with distributed parameters and a computational algorithm for its numerical research are implemented in the form of a computer program. The researches of the model of a polyfrequency vibrating screen with one impactor were carried out. It is shown that a change in the amplitude of the exciting force has a significant effect on the dynamics of the impactor vibrations of a polyfrequency vibrating screen. For such a system, there is an amplitude value at which the impactor passes from the mode without interacting with elastic limiters to the vibro-impact mode. In this case, the movements of the impactor begin to change disproportionately to the change in the amplitude of the exciting force. =56000 N.
3. The researches of a model of a polyfrequency vibrating screen with two impactors were carried out. It is shown that the start of both impactors in a vibrating polyfrequency screen significantly depends on the exciting force, i.e. the problem with the absence of a vibro-impact mode in one or both of the impactors can be solved by increasing the amplitude of the exciting force. In addition, by changing the amplitude of the exciting force, it is possible to achieve chaotic modes of impactors vibration, which leads to oscillations of sieve surfaces with a continuous frequency spectrum, i.e. to modes of the screen operation, which are most appropriate for dehydration and separation of fine fractions of minerals.