CC BY
13
Сучасна практика застосування ei6pa-цшних машин при po6omi з dpiSHUM настням малог ваги зштовхуеться з таким небажаним явищем, як вплив на кнематику вiбpацiйного руху частинок фракцш насгнневог сумШ аеро-динамiчних сил та моментiв. Перюдичний рух повтря вИдносно робочих площин вiбpацiйноi машини виникае за рахунок коливань nакетiв цих площин, як утворюють плост аероди-намiчнi канали. З цього виявляеться акту-альним питания дослИдження процесс взае-модИ робочих оргатв в^рацшног машини з повтряним середовищем з метою обтрунту-вання гх конструктивних доробок. 1снуючи математичт моделi, ят ощнюють параме-три руху повтря вИдносно робочих площин вiбpацiйних машин, дають лише узагальнену картину та являються плоскими. В стат-тi наведено постановку та розрахункову кт-цево-рЬницеву схему розв'язання тpимipног кураевог задачi для обчислення поля швидко-стей та тиску для областi повтря, що зна-ходиться мiж двома паралельними площи-нами, як синхронно коливаються. В задачi використано систему диференцшних piвиянь для опису течп Идеального газу. Ршення кт-цево-ргзницевог схеми зЖйснено за допомогою методу прогонки.
Використання методу прогонки для виpi-шення такого роду завдань дозволяе забез-печити збiжнiсть i стшккть розрахункових схем незалежно вiд кроку i тших паpаметpiв використовуваног сгтки.
Наведено ваpiантрозрахунку, що демон-струе працездаттсть запропонованого методу для заданих крайових умов i параме-тpiв вiбpацiйного режиму роботи машини. Встановлено, що вробочому пpостоpi, укла-деному мiж двома площинами, що коливаються, мають мкце як вертикальна (поперечна), так i горизонтальна (поздовжня) складовi швидкостi руху повтря, як змтю-ються за часом
Ключовi слова: газодинамжа, система диференщальних piвнянь, крайова задача, метод сток, метод прогонки, поле швидко-стей
UDC 519.612:533.69
[DOI: 10.15587/1729-4061.2020.198501|
SOLUTION OF THE SYSTEM OF GAS-DYNAMIC EQUATIONS FOR THE PROCESSES OF INTERACTION OF VIBRATORS
WITH THE AIR
R. Antoshchenkov
Doctor of Technical Sciences, Associate Professor, Head of the Department*
A. Nikiforov Senior Lecturer* E-mail: [email protected] I. Galych Senior Lecturer* V. Tolstolutskyi PhD, Head of Department Department of Automated Control Systems and Engineering Calculations State-Owned Enterprise Kharkiv Machine Building Design Bureau named A. A. Morozov Plekhanivska str., 126, Kharkiv, Ukraine, 61001 V. A n t os h c h e n ko va PhD, Associate Professor Department of Economics and Marketing**
S. Di u ndi k PhD, Associate Professor Department of Armored Vehicles National Academy of National Guard of Ukraine Zakhysnykiv Ukrainy sq., 3, Kharkiv, Ukraine, 61001 *Department of Mechatronics and Machine Parts** **Kharkiv Petro Vasylenko National Technical University of Agriculture Alchevskykh str., 44, Kharkiv, Ukraine, 61002
Received date 19.11.2019 Accepted date 06.03.2020 Published date 24.04.2020
Copyright © 2020, R. Antoshchenkov, A. Nikiforov, I. Galych, V. Tolstolutskyi, V. Antoshchenkova, S. Diundik This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0)
1. 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-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-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.
2. 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 (hydro-dynamic) 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.
3. 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.
4. 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.
pa = pF -grad p,
1
p—2+divV = 0, pc
du + du + du + du 1 dp
dt dx dy dz x p dx
dv dv dv dv 1 dp
--+ u--+ v--+ w— = g---—,
dt dx dy dz p dy
dw ~dt'
dw dw dw
dx dy
-w—= g, -
dz
1 dp p dz '
dp dp dp dp —+u— + v—+ w— dt dx dy dz
du dv dw .
--+ — + — = 0,
dx dy dz
pc
- 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:
u(t)
v(t)
M,B"
M,B"
M,B
= 0,
(1)
(2)
p (t
+p
A = po + V/(t)2(t)2+"/(t)2 (-Sign(t)}),
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; p 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:
p % == po +
+p
v/ (t)2+V/ (t)2+v/ (t)2
sign {V/ (t
sign {V (t)} = \ 1, if V ^
6 1 z W/ 1-1, if Vzk (t)< 0,
(9)
(10)
(11)
(12)
(3)
(4)
(5)
where p is the air density; Vxk (t), Vyk (t), Vzk (t) 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:
I dQ + b dQ+C dQ + D dQ = F ,
dt dx dy dz
(13)
(6)
where
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; gx, gy, gz 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
1 0 0 0 u 0 0 c2
0 1 0 0 , B = 0 u 0 0
0 0 1 0 0 0 u 0
0 0 0 1 0 0 0 u
C =
v 0 0 0 w 0 0 0
0 v 0 c2 0 w 0 0
, D =
0 0 v 0 0 0 w c2
0 0 0 v 0 0 0 w
V (t
C, D, E,G = V K (t),
C, D, E ,G = Po,
(7)
(8)
Q=
u gx
v , f = gy
w gz
pi pc2 0
2
1
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 E, relative to the system of coordinates of the working surface; VK(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; p0 is the atmospheric pressure;
The sweep will run along two axes: the X axis and the Y axis. Thus, for a direct sweep along the Xaxis, 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),
¿=0, ..., b/h, k=0, ..., H/s, is also assigned to the j-th node of the Faxis. 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 F axis are shown by arrows in Fig. 1. Reverse sweep is performed in reverse order.
Z
H/s; k ■
U 1-j
0, 0, 0
0, j, H/s
b/h,j,0
X
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, t=1,..., T, one computes the elements of Xj,T and j tensors, using the following recurrent ratios:
X
(j+1),T "
j+1),(T-1)/
-AA°,(t-1) (AAA(+
+ AAA( j-1) (T-1) ( AAA(+
j+1),(T-1) j X j,T
y( j+1),T - - AAA(-;-1)(T-1) ( AAA(+j+1),(T-1) ) X
XX(j+1),t,t - (t-1) ( AAA(+j+1),(T-1) j X
where
(14)
(15)
(j+1),T'
AA0T-
jT
-Aj A0 j t A+, j ,t 0 0 0 -A-j,T A0,j,t a+,j,t 0
0 0 0 0 0
0 0
A+.
, (16)
0 1/2h - b^
_ j,T _ 1/2h - Bi,; ,1,T N 1/2h - bi,;,,t
1/2h - Bi,;№,t 0
A° j ,T -
-1
"2S Di, j ^ 0
1 1
A1 2s D,j^
-1
~2s D,j"T-1)
J_
At
I
2 s
i, j,3,T
(17)
0 ... 0 ...
Di,j,k,% -
ui, j,k,T 0 0 c2
0 ui, jkr 0 0
0 0 ui, jkT 0
0 0 0 ui, jkr
®i,j,k,T 0 0 0
0 Wi, j,k,T 0 0
0 0 Wi, j,k,T c2
0 0 0 w- i, ],k ,T
AAA+T-
AA+
0
1/2l -Ci,;
1/2l -C
'i, j ,Hjs,T
C
, j,k ,T 0 0 0
0 vi, j,k,T 0 c2
0 0 vi, j,k,T 0
0 0 0 vi, j,k,T
Fj,T-
F1,j
bjh, j;
F -
i, j ,T
Fr + At , j 0, 1
Fr + At IQ<,j'H'sr
Ft -
gx (rAt) gy (rAt) & (rAt) 0
(19)
(20)
0 aa+j,r
aa+jr , aaa-t - aa+jr , (21)
AAb/ h, j ,t 0
(22)
(23)
(24)
The computed elements are saved for all the steps j=0, ..., a/l.
The elements of the tensor-section for 0 the far-right position (at the end of the di-
rect sweep), Sa/;T are determined:
0 ¥„.,
0 0 ...
f Dh
2s ¿.I —
, (18)
T
T
b I ha! l,r
v
h
0
0
0
where
ia¡l,T
Qia/l,0,T
Qi,a/l,l, T
i¡l.HjsT
q.,
i/l,k,T "
@¡i,k,T
,a¡l,k,T
,a/l,k,T
Pia/lk,*! pC
t=1/2fl) (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 inhibited - the boundary condition (9). Its relative velocity becomes zero. The kinetic energy of air movement passes into the energy of excess pressure - the boundary conditions (10) to (12).
The values for m^//^, via/i,k¡l, wia/i,k¡l, pi,a//,k,z 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:
-(j-1),T = + yj,T,
j = a//, (a//-1),(a//-2),... ,1. (25)
The field of velocities and pressures is computed.
5. 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 i=0 (Fig. 2), í=1/4H (Fig. 3) and £=1/2H (Fig. 4). The time point i=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 í=1/4 H 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 í=1/2 H are identical to the time point i=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 i=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
Fig. 2. The field of velocities and pressures for time point t=0
Fig. 3. The field of velocities and pressures for time point t=1/4H
Fig. 4. The field of velocities a nd pressures for time point
i=1/2n
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
u
v
w
reached near the lower plane. When moving from bottom to top: near the top plane.
6. 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/20], 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 i=1/2fi (Fig. 4) mirrors the pattern for t=0 (Fig. 2). It is obvious that over the full period of oscillations, for i=1/fl, 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.
7. 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).
References
1. Adamchuk, V. V., Prilutsky, A. N. (2017). Theoretical studies of non-grid pneumovibrating centrifugal separation of seed mixtures. Mekhanizatsiya ta elektryfikatsiya silskoho hospodarstva, 5 (104), 17-27.
2. Khodabakhshian, R., Emadi, B., Khojastehpour, M., Golzarian, M. R. (2018). Aerodynamic separation and cleaning of pomegranate arils from rind and white segments (locular septa). Journal of the Saudi Society of Agricultural Sciences, 17 (1), 61-68. doi: https://doi.org/10.1016/j.jssas.2016.01.003
3. Posner, E. S. (2009). CHAPTER 5: Wheat Flour Milling. WHEAT: Chemistry and Technology, 119-152. doi: https:// doi.org/10.1094/9781891127557.005
4. Nazir, A., Khan, K., Maan, A., Zia, R., Giorno, L., Schro n, K. (2019). Membrane separation technology for the recovery of nu-traceuticals from food industrial streams. Trends in Food Science & Technology, 86, 426-438. doi: https://doi.org/10.1016/ j.tifs.2019.02.049
5. Okunola, A. A., Igbeka, J. C. (2009). Development of a reciprocating sieve and air blast cereal cleaner. African Crop Science Conference Procedding, 9, 3-8.
6. Luk'yanenko, V. M., Galich, I. V., Nikiforov, A. A. (2015). Mehatronnaya vibratsionnaya semyaochistitel'naya mashina. Mehani-zatsiya sel'skohozyaystvennogo proizvodstva, 156, 413-419.
7. Dzyuba, O., Dzyuba, A., Polyakov, A., Volokh, V., Antoshchenkov, R., Mykhailov, A. (2019). Studying the influence of structural-mode parameters on energy efficiency of the plough PLN-3-35. Eastern-European Journal of Enterprise Technologies, 3 (1 (99)), 55-65. doi: https://doi.org/10.15587/1729-4061.2019.169903
8. Antoshchenkov, R. V. (2017). Dynamika ta enerhetyka rukhu bahatoelementnykh mashynno-traktornykh ahrehativ. Kharkiv: KhNTUSH, 244.
9. Chattot, J. J., Hafez, M. M. (2015). Theoretical and Applied Aerodynamics. Springer. doi: https://doi.org/10.1007/978-94-017-9825-9
10. Sirca, S., Horvat, M. (2018). Computational Methods in Physics. Graduate Texts in Physics. Springer. doi: https://doi.org/ 10.1007/978-3-319-78619-3
11. Strauch, D. (2009). Classical Mechanics. Springer, 405. doi: https://doi.org/10.1007/978-3-540-73616-5
12. Volkov, K. N. (2009). Preconditioning of the Euler and Navier-Stokes equations in low-velocity flow simulation on unstructured grids. Computational Mathematics and Mathematical Physics, 49 (10), 1789-1804. doi: https://doi.org/10.1134/ s0965542509100133
13. Durst, F. (2008). Fluid Mechanics. Springer, 723. doi: https://doi.org/10.1007/978-3-540-71343-2
14. Dimov, I., Fidanova, S., Lirkov, I. (2015). Numerical Methods and Applications. Springer, 313. doi: https://doi.org/10.1007/978-3-319-15585-2
15. Mastorakis, N., Sakellaris, J. (Eds.) (2009). Advances in Numerical Methods. Lecture Notes in Electrical Engineering. Springer. doi: https://doi.org/10.1007/978-0-387-76483-2
16. Lei, L., Chaolu, T. (2011). A new method for solving boundary value problems for partial differential equations. Computers & Mathematics with Applications, 61 (8), 2164-2167. doi: https://doi.org/10.1016/j.camwa.2010.09.002
17. Zavgorodniy, A. I., Sinyaeva, O. V. (2012). Dvizheniya shara v vozdushnom potoke mezhdu vibriruyushchimi ploskostyami. Vi-bratsiyi v tekhnytsi ta tekhnolohiyakh, 3, 20-27.
18. Lukianenko, V., Nikiforov, A., Petrik, A. (2017). The method of calculating the velocity field of the air environment between two equidistant ploskostyami when they commit synchronous harmonic. Vestnik HNTUSG im. P. Vasilenko, 180, 100-105.
19. Luk'yanenko, V. M., Nikiforov, A. A. (2017). Postanovka zadachi rascheta polya skorostey vozdushnogo sredy mezhdu dvumya ekvidistantnyh ploskostyami pri osushchestvlenii imi sinhronnyh garmonicheskih kolebaniy. Inzheneriya pryrodokorystuvannia, 2, 33-37.
20. Chadwick, P. (1999). Continuum Mechanics: Concise Theory and Problems. Dover Books on Physics, 192.
