Solution of the System of Gas-Dynamic Equations for the Processes of Interaction of Vibrators With the Air

The modern practice of using vibratory machines involving small seeds of low weight faces such an undesirable phenomenon as the effect exerted on the kinematics of vibrational movement of particles of fractions of the seed mixture by the aerodynamic forces and momenta. The periodic movement of air relative to the working planes of a vibratory machine arises due to fluctuations in the packets of these planes, which form flat aerodynamic channels. Consequently, the issues of studying the processes of interaction between the working bodies of vibratory machines and the air environment, aimed to justify their structural improvements, appear relevant. Existing mathematical models, which assess the parameters of air movement relative to the working planes of vibratory machines, produce only a generalized pattern and are flat. This paper proposes a statement, as well as an estimated finite difference scheme, of solving a three-dimensional boundary value problem on calculating the field of velocities and pressures in the region of air, located between two parallel synchronously oscillating planes. The problem employs a system of differential equations to describe the flow of the perfect gas. The finite difference scheme has been solved by a sweep method.<br><br>Using the sweep method to solve these kinds of problems makes it possible to ensure the convergence and stability of estimation schemes, regardless of the step and other parameters of the grid applied. A variant of the calculation has been given, which demonstrated the feasibility of the proposed method for the assigned boundary conditions and parameters of the vibrational mode of machine operation. It has been established that in the working space enclosed between two oscillating planes there are both vertical (transverse) and horizontal (longitudinal) components of air velocity, which change over time.


Introduction
The aerodynamic properties of seeds are widely used in seed cleaning practice, especially for the removal of impurities (straw, spike particles, etc.). Such impurities are sharply different in terms of specific weight from seeds; they have different critical velocity and are easily detached in the airflow. This principle of cleaning is fully and successfully used in simple (winnows) and complex seed-cleaning machines. Critical seed velocity depends to a large extent on their shape: in a spherical seed, it is quite constant, so it can be employed during cleaning. There is a group of weed seeds that vary in their sailing ratio; they can be completely removed by exploiting this property.
The effective separation of seed mixtures with pronounced aerodynamic properties is carried out using devices that separate them based not on a single attribute, but rather based on a set of physical-mechanical properties: shape, roughness, and elasticity [1,2]. These devices include vibratory machines, which have shown high separation efficiency for many small-seeded crops [3][4][5].
For such machines, it is of interest to investigate the process of interaction between the working surfaces of the machine and air. The need to use gas-dynamic models is caused by that there is an air movement in the airspace between the oscillating work surfaces. This phenomenon has a significant impact on the pattern of vibrational movement of seeds with pronounced aerodynamic properties. To investigate this interaction process, it is necessary to have adequate mathematical models for predicting airflow parameters over vibrating work surfaces, depending on the parameters of the operational mode and design features of vibratory machines.
The current theory of the vibrational movement of small seeds (particles) does not fully take into consideration the aerodynamic factor. This is mainly due to the lack of three-dimensional air movement patterns under the influence of the working organs of a vibratory machine.
This inhibits further improvement of vibrational cleaning methods and means as the most effective way to separate small-seeded crops. For example, a new mechatronic vibratory cleaning machine [6][7][8], while providing significant performance improvement, needs to be refined to compensate for the aerodynamic factor. The selection of rational structural parameters of the improved machine is possible only based on multivariate studies on modeling the workflow taking into consideration the dynamics of air mass movement. The kinematic parameters of air medium movement that interacts with the processed seed crop must be calculated using gas-dynamic equations that could be solved by applying modern numerical methods. It is desirable that the method to be used, while producing a three-dimensional pattern of air movement, should at the same time not require an excessive increase in computational resources, as is the case, for example, in the gas-dynamic calculations of thermal machines.

Literature review and problem statement
The numerical methods for solving the problems of hydro-gas-dynamics are constantly evolving and are applied in many practical areas of activity. These areas are related to the design of aircraft, water, and ground-based vehicles, various assemblies and devices, whose operation implies taking into consideration the impact of the surrounding air (gas) or water environment. Up to now, a large number of estimation schemes and models have been constructed, allowing the calculation of parameters of the gas-air (water) environment when interacting with the structural elements of designed vehicles (assemblies) [9]. The applied estimation approaches are mainly based on the grid method, the method of generators, or the Massot method [10,11]. These methods, while demonstrating undeniable advantages in terms of simplicity and versatility, make it possible to resolve the issue of the non-linearity of differential equation systems, which typically describe the examined gas dynamic (hydrodynamic) processes. The price incurred is the instability and unsatisfactory convergence of the solutions derived, which depend on the technique and grid parameters for splitting the regions under study.
The problem of convergence of the numerical solution to gas-dynamic equations is central to the research of estimation schemes used in the field of thermo-gas-dynamic processes.
An iterative method for solving the Euler's finite difference equations was proposed in [12]. Given the simplicity of the method, the solution process has a slow convergence rate, which is unacceptable for variable calculations.
In [13], to improve convergence, a matrix time step, and a method of directional coarsening of the grid were proposed. That made it possible to significantly increase the speed of convergence of calculations but significantly complicated the algorithm.
Papers [14,15] examine the methods of multi-level multiple grids, where the error of the solution obtained on a small grid is transferred to a large grid, and then the smoothed solution, obtained on a large grid, is transferred back to a small grid. The methods proposed, while outperforming the previous method in terms of convergence rate and the achieved accuracy of calculations, are even more complex and costly in terms of the consumption of computational resources.
It was noted in [11] that for the case of a perfect gas when differential equations are brought to a quasi-linear form, thereby reducing the estimation model to a boundary value problem, the toolset of applied computational methods could be supplemented with a sweep method. According to the studies reported in [16], it is claimed that a given method is insensitive to how a region splitting grid is formed. When decreasing a breakdown step, the accuracy of the solution always improves, which converges to a certain value. To ensure the stability of the solution, no additional measures are required, for example, the use of data formats with an increased number of bits.
The statement of a calculation scheme for solving the boundary value problem by a sweep method for a two-dimensional case when solving a problem on heat exchange on the plane was proposed in [9]. As regards the case of a three-dimensional system of gas dynamics equations, the authors are not aware of any results in obtaining the estimation schemes to solve them using a sweep method.
For the case of studying relatively simple gas-dynamic processes: the use of the perfect gas model, within an acoustic range of velocities and pressures, it is advisable to consider the sweep method as a priority. This is due to a significant gain in the rate of computing at a relatively simple analytical development of the computational scheme.
Thus, when considering the processes of interaction between the working surfaces of vibratory machines and an air environment, when there are weak (acoustic) disturbances, it seems appropriate to conduct research on the analytical development of estimation schemes using a sweep method.

The aim and objectives of the study
The aim of this study is to devise a method for solving a boundary value problem for a three-dimensional system of differential equations of air dynamics under the influence of the working bodies of a vibratory machine using the sweep method. Applying the sweep method to solve these kinds of problems makes it possible to ensure the convergence and stability of calculation schemes, regardless of the step and other parameters of the grid used.
To accomplish the aim, the following tasks have been set: -to construct a numerical sweep algorithm for three orthogonal axes; -to implement the built computer algorithm on PC using the MATLAB programming environment, in order to demonstrate its feasibility by computing the field of velocities and pressures of air medium for the characteristic positions of the working bodies of a vibratory machine.

Construction of a numerical sweep algorithm for three orthogonal axes
In a general form, the statement of a boundary value problem to calculate the field of velocities and air pressures, located between two parallel synchronously oscillating working planes of a vibratory machine, is given in [19]. The results are presented in the analytical and finite difference form.
The Euler equation, supplemented with a continuity equation, was used as a mathematical model of the process under study.
where a is the air medium acceleration vector; F is the vector of acceleration due to the action of mass forces (gravity); p is the air pressure at the point in question; ρ is the air density, V, p  is the motion velocity vector and the velocity of air medium pressure change at the point in question, respectively; c is the speed of sound.
In the coordinate form, the system of equations being solved takes the form: 1 , where u, v, w are the projections of the air medium velocity vector, V, onto the X, Y, Z axes of the selected coordinate system; g x , g y , g z are the projections of the acceleration of free fall onto the axes of the selected coordinate system. The boundary conditions for the borders of the region confined between two synchronously oscillating work surfaces can be recorded as follows: -for edges: C, D, E and G (along the contour of the calculated area) there is an unimpeded relative movement of air and the pressure is equal to the atmospheric pressure 0 , , , , where V(t)/C, D, E, G is the velocity vector of air particle motion, which belong to the C, D, E and G edges of region Ξ, relative to the system of coordinates of the working surface; V K (t) is the velocity vector of the oscillations of points at the working surface relative to the inertial system of coordinates; p/C, D, E, G is the air pressure along the C, D, E, and G border; p 0 is the atmospheric pressure; -for the A and B edges, which come into contact with the surfaces of the lower and upper working planes of a vibratory machine, there is complete braking of the air when it comes into contact with the surface (the air at rest is set into motion by the oscillating working surfaces). A positive or negative pressure difference Δp is formed. The sign of this difference is determined depending on the motion direction of the working surface. That is, the boundary conditions for edges A and B take the form: where ρ is the air density; are the projections of the oscillation velocity.
An estimation scheme was proposed in [20], which implements the sweep for three orthogonal axes of the coordinate system associated with the working planes of a vibratory machine.
To derive an estimation scheme, the system of equations (3) to (6) is written in a matrix form: The sweep will run along two axes: the X axis and the Y axis. Thus, for a direct sweep along the X axis, each j-th node belonging to the Y axis is assigned with a set of nodes lying on the vertical axis, which passes through the j-th node. The formed left boundary of the studied region along the j-th section is swept to the right boundary by moving it from the ZOY plane along the OX axis. The resulting set of nodes (i, j, k), i=0, …, b/h, k=0, …, H/s, is also assigned to the j-th node of the Y axis. Next, the formed section is swept to node j=a/l (to the side end of the studied region). The directions of direct sweep along the X axis and Y axis are shown by arrows in Fig. 1. Reverse sweep is performed in reverse order.

Fig. 1. Sweep scheme
Papers [19,20] give the finite difference notation of the system of equations (16) and its transformation aimed to derive the recurrent direct-to-reverse sweep ratios.
In the course of a direct sweep, for j=0,…, a/l, τ=1,…, T, one computes the elements of X j,τ and y j,τ tensors, using the following recurrent ratios: ( 1), 1    The values for u i,a/l,k,τ , v i,a/l,k,τ , w i,a/l,k,τ , p i,a/l,k,τ are determined based on the boundary conditions for edge D.
Next, we run a reverse sweep. We compute the tensor-section elements for all intermediate positions to the left of the right end of the sweep interval: The field of velocities and pressures is computed.

The results obtained and their analysis
Based on the proposed analytical expressions (14) to (25), we have constructed an estimation algorithm, which is implemented in the applied software package MATLAB designed to solve the problems of technical calculations. The results that were obtained using it are given for time moments t=0 (Fig. 2), t=1/4Ω (Fig. 3) and t=1/2Ω (Fig. 4). The time point t=0 corresponds to the neutral position of a vibratory machine's planes. There is a maximum (by module) value of the velocity magnitude and a zero value of the acceleration of the movement of oscillating planes. The time point t=1/4Ω corresponds to such position of planes where their deviation from the zero position is maximum. There are a maximum acceleration and a zero velocity of plane movement. The motion parameters and the position of the planes corresponding to point time t=1/2Ω are identical to the time point t=0, but the movement speed of the working planes of a vibratory machine here is directed in the opposite direction.
The illustrations above (Fig. 2-4) show that in the working space confined between two oscillating planes there are both vertical (transverse) and horizontal (longitudinal) components of air velocity, which change over time. The law of change in the longitudinal and transverse components of velocity is periodic, with a change in the movement direction. The superposition of these two movements produces a complex distribution pattern of the movement speed of the elements within the studied air continuum, where there are uneven velocities both vertically and horizontally in the estimation area.
The resulting distribution of pressures is characterized by unevenness for the height of the region. When moving from top to bottom, at time point t=0 (Fig. 2) there is a maximum pressure difference at the inner surface of the upper plane. In an extreme position of the working bodies, when the velocity of planes is zero, there is no dynamic pressure difference (Fig. 3). When moving backward (time point t=1/2Ω) (Fig. 4), there is excess pressure at the inner surface of the lower plane.
The magnitude of excess pressure corresponds to the magnitude of relative air velocity. The layer of air directly adjacent to the inner surface of the plane is inhibitedthe boundary condition (9). Its relative velocity becomes zero. The kinetic energy of air movement passes into the energy of excess pressurethe boundary conditions (10) to (12).  Thus, the maximum air velocity relative to working surfaces is achieved in a position where the working bodies of a vibratory machine pass the neutral position. When moving from top to bottom: the maximum relative air velocity is reached near the lower plane. When moving from bottom to top: near the top plane.

Discussion of results of studying the application of a sweep method
The above calculation results (Fig. 2-4) demonstrate that the field of velocities and pressures changes in accordance with the harmonic law, over a period equal to the period of oscillations of the working planes of a vibratory machine. For t=[0; 1/2Ω], which is equal to the half-period of oscillations of the working planes, the direction of air motion changes to the opposite, and the pattern of the distribution of velocities and pressures for t=1/2Ω (Fig. 4) mirrors the pattern for t=0 (Fig. 2). It is obvious that over the full period of oscillations, for t=1/Ω, the pattern of the distribution of velocities and pressures would take the form shown in Fig. 2 for the starting time point t=0. The movement of the air occurs under the influence of pressure drop caused by the dynamic pressure, which is exerted on the air by the incoming plane, on the one hand, and by sucking near the outgoing plane, on the other hand.
The resulting pattern of the dynamics of air movement does not include the vortex phenomena that are likely to actually occur. However, the model used, which does not take into consideration the viscosity of the air, does not make it possible to derive swirls in a calculated way. When using, as a kinematic model, instead of the Euler equation, the Navier-Stokes equation, such vortex effects would be obtained. However, at the same time, the estimation scheme would be much more complicated and, perhaps, would not make it possible to implement the sweep method.
However, for the practical aspect of this issue, the noted methodical flaw is of inconsequential importance. To study the effect of the aerodynamic factor on the character of the vibrational movement of seeds, the most significant is the accounting of the tangential components of air velocity. And the proposed calculation method makes it possible to compute them successfully, and at the low time and machine memory costs.
Further prospects for the development of the results reported here imply the improvement of the formalized description of the boundary conditions, which make it possible to take into consideration the various structural elements used to eliminate the harmful effects of the aerodynamic factor. In addition, we believe, the proposed algorithm has an independent value for its application in the field of gas dynamic calculations over the acoustic range of gas currents. The convergence of the solution automatically provided by the sweep method could significantly simplify the various computational problems on determining the parameters of gas (air) in subsonic flows, without taking into account the viscosity.

Conclusions
1. An algorithm for sweeping along three orthogonal axes has been constructed in order to solve a system of gas-dynamic equations recorded for the case of perfect gas over the acoustic range.
2. The algorithm calculates the field of velocities and pressures in the air mass, located between two parallel synchronously oscillating planes of a vibratory machine. The resulting three-dimensional air dynamics pattern does not require significant computational resources, so it is advisable to use it to build complex models to study the processes of interaction between the air environment and the oscillating elements of vibratory machines when solving the tasks of their design (improvement).