MSC 65N06, 76G25
DOI: 10.14529/mmp170411
MATHEMATICAL AND SOFTWARE SUPPORT
FOR 3D MATHEMATICAL MODELLING OF THE AIRFLOW IMPACT ON THE OPTICAL-MECHANICAL UNIT MOUNTED IN THE AIRCRAFT UNPRESSURIZED COMPARTMENT
I.E. Ivanov, I.A. Kryukov, E.V. Larina, V.L. Miroshkin
Moscow Aviation Institute (National Research University), Moscow, Russian Federation E-mail: [email protected], [email protected], [email protected], [email protected]
The problem of three-dimensional mathematical modelling of the effect of air flow on an optical-mechanical unit (OMU) located in the unpressurised compartment of the aircraft, is considered. To solve this problem, a mathematical model of gas dynamics based on the solution of a complete system of Navier - Stokes equations that describe the dynamics of a turbulent, spatially unsteady flow of a viscous gas is developed. The software for simulating the process of flow past a WMU model in the aircraft compartment was created. The effect of the air flow on the OMU is described by the torque acting on the OMU from the airflow side. A numerical method for solving the three-dimensional gasdynamic problem is presented. The numerical method is based on the numerical high order Godunov scheme, realized on an irregular grid with arbitrary cells (tetrahedral, prismatic shape). Flows of conservative variables are calculated by solving the Riemann problem with an approximate AUSM method. The system of equations is supplemented by a two-parameter k-model of turbulence, modified for the calculation of high-speed compressible flows. To significantly reduce the cost of computing resources, it is suggested to use stochastic models of the effect of air flow on WMU. A general simulation algorithm is described.
Keywords: optical-mechanical unit; mathematical modelling; the effects of air flow.
Introduction
For a number of tasks on tracking of space objects in optical wavelength range it is possible to use optical-mechanical unit (OMU), mounted in the aircraft's fuselage, the OMU can be placed only inside the fuselage on the special platform. The fuselage has a special window through which the OMU performs tracking of the space objects. The output window is sealed against a direct flow by the special radome. Due to the flowing of the airflow through the niche the OMU fairing arise turbulent flows, leading to pressure fluctuations at the place of installation of OMUs. Pressure fluctuations cause vibrations that may degrade the accuracy of the OMU's work. Before the field experiments (flight testing) the only way to assess the impact of airflow on the OMU is a mathematical modelling. The mathematical modelling of the airflow allows to simulate changes in values of the airflow major parameters in any prescribed point of the OMU. However, the numerical solution of gas dynamics equations requires enormous computational resources that makes it almost impossible to express analysis of various design and technological solutions [1]. Therefore, a gas dynamic model is appropriate to use as a stage of data preparation for the subsequent identification of the stochastic model [2]. In this article the advanced approach to the mathematical modelling of the impact of airflow on the OMU is proposed.
1. Mathematical Gas-Dynamics Model
The appearance of unsteady flows in cavities and compartments containing discrete components in the spectrum of fluctuations is associated with the development of large-scale structures in the mixing layer, and the feedback from the source of vibrations and place of formation of large-scale vortices. The formation of ordered oscillations at discrete frequencies may contribute to the resulting pressure wave transmitted upstream.
Thus, the physical model of the flow should be a model of a viscous compressible gas with the formation of large-scale structures that interact with the acoustic field and shock waves (if any), propagating against the flow of gas.
This physical model is adequately represented by the mathematical model based on the system of multi-dimensional non-stationary Xavier Stokes equations for compressible gas with appropriate initial and boundary conditions and closing relations.
The complete system of Favre averaged Xavier Stokes equations, describing the spatial unsteady compressible flow, viscous gas and used in mathematical modelling, refer-to [1,3 7].
The geometric model of the OMU is presented in Fig. 1.
Fig. 1. The geometric model of the OMU
2. Numerical Method for 3D Gas Dynamics Modelling
To solve the equations the volume control method, namely quasimonotone numerical scheme, is used. This method is a modification of the circuit S.Iv. Godunov [8 10]. Volume control volume are based on application of the Gauss - Ostrogradsky theorem. According to this theorem for the vector field $ = ( ) and a closed volume V, wich is
bounded by the outer surface S with the normal vector vecn, the equality is fulfilled
JJJ div ($^jdV = j) $ • ndS, (1)
V S
where the divergence of the vector field in the Cartesian coordinate system have div ^ $ >dXL- According to (1), for the system of equations [1] and each index number corresponding
to the equation system integrals over the volume of the control of non-viscous and viscous streams may be converted:
fff( dx+f+%) -=/ =l. «
V S
jam+g+h) dv=/ * ■-=m. (3,
V S
where the components of the vectors of vicous flows F^ Gj Hj, standing consistently
—T
in every j-th equation systems are joined into the vector Lj = ( Fj Gj Hj ) , and the components of the vectors of viscous flows Fvj, Gvj, Hvj, facing successively in each J-th equation systems are joined into a vector Mj = ( Fvj Gvj Hvj )T. Thus, change of the vector of conservative variables in each cell is defined by flows across the border, and by sources / sinks in volume, which are defined by volume integrals from source component members in vector. The numerical scheme is reduced to the assumption that the distribution of variables within a cell, finding flows through the faces of the cells and the integration time.
The differential-difference representation of the original system equations written for the calculation of the cell is given by:
d(q(t) + R (t)) 1
dt V
Y.L f-i + Mf-i . fi
(4)
where V is cell volume, Lfi is flux through the face Fi, corresponding to the formula (2), MMfi is flux through the face of fi, according to the formula (3), the summation carried out on all limiting the cell faces;
q(t) = v J q(x,y, z,t)dxdydz (5)
is average over the entire volume of the cell value of the vector of conservative variables
qt
similar formula R (t).
Fig. 2 presents the types of cells that are valid for conducting numerical experiments. Splitting the computational domain into control volumes should be carried out in such a way that adjacent cells together had only one common face. Thus, the number of adjacent cells is determined by the number of cell faces. Cells having a common face each other, are called adjacent.
Splitting the computational domain into control volumes should be carried out in such a way that adjacent cells together have only one common face. Thus, the number of adjacent cells is determined by the number of cell faces, and cells having a common face are called neighbors. Non-viscous flows on each face of the control volume (cell) corresponding to the equations of conservation of mass, momentum and energy, are calculated using the approximate solution of the Eiemann problem, namely by the method AUSM-M (AUSM -advection upstream splitting method [4]). AUSM-methods for calculating viscous flows are quite economical and suitable for calculating viscous flows. At this move to a coordinate
Fig. 2. Types of acceptable computational cells (control volumes): a) the pyramid ABCDE, which lies at the base of arbitrary the quadrilateral; b) a prism ABCDFE formed as a result of pulling (parallel transport) along the triangle ABE a side edge, such as along BC; c) the tetrahedra ABC; d) ABCDEFHG hexagons which edges are arbitrary the quadrangles
system is associated with the face and we consider only the part of velocity which is perpendicular to the face, sensing tangential velocity as a portable scalar. We get a vector of flow
/
F
\
pu \ / paM \ f 0 \
puu + p puaM V
pvu = pvaM + 0
pwu pwaM 0
phu / \ phaM \0)
(6)
where h = ^ is gas enthalpy, u is the projection of velocity on the direction normal to the face (u = un). Valu es p, u, v, w, h are are considered as passive scalars transported by convective speed M and rely equal to the corresponding values from the "left" cell if the convective velocity on a face directed toward the "right" cell and vice versa. Value p represents the action of the pressure.
F
is the sum of two components. M is the Mach number obtained in AUSM-M method as the sum of the "convective velocity", partly coming from the "left" (M+) and partly from "right" (M— cell (M = M+ + M—). Part of the flow vector, associated with pressure, is also decomposed into the sum of two components from the "left" (P+) and "right" (P-) cell in dependence of the Mach numbers and pressure in these cells (p = P+ + P—-The splitting of the Mach number and pressure is carried out using the following:
Mt/R
(M )
{
± 4 (m ± i)2
2 (M ±\M |)
|M| < 1, \M\ > 1,
(7)
|M| < 1, |M| > 1.
(8)
Under the assumption of a constant distribution of the parameters within the cells method is only the first order accuracy in space. To achieve the second order accuracy is used piecewise-linear recovery, that is, within a computational cell is assumed linear distribution of the physical parameters. In this case, vectors of variables on the left and right faces of the cell, which separates adjacent cells with numbers textit i and textit j, necessary for solve the problem of decay of a discontinuity can be defined as follows
where q is a scalar variable, Vq - gradient and the variable r is a vector extending from the center of the cell in the center of the face.
It is well-known that the recovery of the second or higher order requires the use of limiters for suppressing artificial oscillations of solutions in the fields of high gradients. We use Barth and Jespersen limiter, specially designed for unstructured meshes. The value of the limiter in each cell is selected as a minimum limiter value determined for each face which bounds the cell. The value of the limiter for the face in this embodiment is determined by the condition that the values on the face do not exceed the maximum and minimum values in the cells of the whole pattern, consisting of counterflow cell and all its adjacent cells. The resulting limiter is of the form:
where function $ (x) = min(x, 1), index f corresponds to the value on the face of the cell, P corresponds to the opposite cell, N corresponds to neighboring cells of P, ^^corresponds to the radius vector of the center faces, rP corresponds to the radius vector of cell center PP
After determining the value of the limiter for the cell, gradients, required for the linear recovery in the calculation of flux through the faces, are multiplying by the limiter value found for the corresponding cells.
To calculate gradients the least squares method is used. From template, which includes some current cell C and adjacent to her cell nb, is held hyperplane that restores distribution of a quantity Q within the cell with the maximum order of accuracy:
qL = qi + vqi • rL qn = qj + ^qj • tr '
(9)
(10)
фр = min (фf ).
Enb Axlb J2nb Axnb&Vnb Ysnb ^xnb^znb
J2nb AXnb^Vnb J2nb AVnb J2nb AVnb&znb
J2nb AXnb^Znb J2nb ^Vnb&Znb J2nb Azrt>
'nb
where Axnb = xnh - xc, Aynb = Vnb - Vc, Aznb = Znb - zc, Aqnb = qnb - qc, the slope of the plane corresponds to the required values gradients within the cell (qX° = dx |c,
qC) = dq
qy = dy
q(c) = §q | X
c' qz = dz\C>
The gradients of velocity and temperature, and turbulent values on the faces of the cells necessary to calculate the viscous flows may be calculated as the average of the calculated gradients of the centers of the cells:
Vqf • n = 1 (VqL + Vqu) • n.
However, it is known that such an approach can lead to a mismatch decisions on rectangular or hexagonal grids, so we use a modification of the formula [5]:
Vqf • n =7^-aLR + 1 (VqL + Vqu) • ( n - aLRs) ,
Wtr - r l\\ 2
n s
cells, WrR - rL\\ is the distance between the centers of the cells L and R, aLR is the scalar production aLR = s • n, index f corresponds to the value on the cell's face.
For the time discretization we use the explicit Eunge - Kutta method of the second-order accuracy:
iil) = Y,almQT+PimAtLm], i = i,2,...,p, (12)
m=0
which flows in the cell c are the sum of the convective and diffusive fluxes through the faces of the cell LCm) = L(qt^ | ) + M^^ Aj ), while the indices
of adj1, ... , adjn correspond to adjacent cells to the cell c (i.e., having a common face), the indexes l and m correspond to the temporary layer. The order of accuracy, as well as
preservation TVD properties achieved by selecting an appropriate set coefficients alm, ¡3lm p
tn
point in time tn+1 = tn + At. We use Eunge - Kutta method of the second order with the following set of parameters: p = 2, a10 = 1, (310 = 1, a2o = 1/2, a21 = 1/2, = 0, ¡321 = 1/2. Time step At is calculated taking into account the non-viscous and viscous restriction on the step size. Approximation of the system of the last two levels is carried out by the control volume method:
HI it (pk) - (Pk - pe + pD)dV = (d(-pk) - (Pk - pe + pD)) fff dV = vc vc
= If (-£ (puk) - dy (pvk) - d (pwk) + £((p + fk) fx)+ (13)
+l((P + f,) dt) + £ ((P + t) £ ))dS,
where the second inequality is written by the mean value theorem. Then the equation has the form
Vol • pki3+ ~t pki = -CONyn + DIFn + Vol • PZ+1 - Vol • pej1,
where Vol is the volume of cell VC CONVn is the part of the equation responsible for explicit approximation of convective members, and DIFn is the part of the equation
responsible for the explicit approximation of the diffusion members. The production of kinetic energy can be approximated explicitly P'n+1 ~ Pkij>, since the member dissipation of kinetic energy of the turbulence presents in the form of
-Vol • At • pj _ -Vol • At • pj j « {-Vol • At • pj kn+l.
Thus, we have the equation for calculating kinetic energy turbulence in a new temporary layer
m , л 4 (m , DIFli-CONVn\
n+l _ Pknj + At ^ (Pkij + Voir j
j j + At • p Щ
г3
Similarly, the equation for dissipation is approximated as follows:
n -L Л+ (r fn pn n I DIFn-CONVij \ n+l P£ij + At • \T£ C^lflijPkij kn + Vol J
£
ij
j + At • Tl C£2f2ijp щ
3. Assessment of the Airflow Impact on the OMU
The impact of time-varying air pressure torque on the OMU can lead to strong vibrations and even resonance oscillations of the OMU, that is totally unacceptable. Therefore, the time pressure forces defines the dynamic loading on the the OMU. The time pressure forces are calculated relative to the unit geometric center (Fig. 1).
Forces and torque are determined by the formulas:
N N
Fx Fix PiASiUix,
i=l i=l NN
Fy Fiy piASiniy J
i=1 i=1 NN
Fz = J2 Fiz = J2 PiASiUiz,
i=1 i=1 N
MM = E ri X F,
i=1
where the summation is over all the wall cells, pi is the pressure in the parietal cell, ASi is a surface element (along the solid surface), on which the force acts, щ = ( uix Niy uiz ) is a unit normal vector to the surfaces. Because of the symmetry Fz = 0. щ is the radius vector drawn from the center of the the OMU to the area middle of the unit. is local
Fi = Fix Fiy Fiz only component of the torque about the axis Z.
Because of the symmetry with respect to the xy plane of the coordinate system (Fig. 1), corresponding components of the air pressure torque are equal to 0. The change in time Z
The time interval corresponds to approximately 0,08 s time transients in the calculation scheme. Dynamic loading on the OMU is negligible, of the order ±60N • m. However, even such a small stress can result in significant impacts on the unit in the event of resonance. The impact of pressure forces on the OMU should be simulated to analyze the drive dynamics on the possibility of resonance.
Fig. 3. The torque about the OMU center. MZ is the component of the torque along with the z axis direction, time is in seconds (x-axis), toque is in N • m (y-axis)
4. Description the Simulation Algorithm
Theoretically, the proposed gas-dynamic model in the presence of the required computer capacity allows you to simulate the impact of airflow on the OMU with any given accuracy. But when you run a simulation on the proposed gas dynamic model, the computation time for the simulation was about 24 h per one period of oscillation at a frequency of 50 Hz. Using time series, it was found that for the end of the transition process in at simulation requires at least 0,1 s. It is clear that to carry out the simulation of gas dynamics for a few seconds is almost impossible. It is therefore proposed to use the short time interval to simulate gas dynamics and to use its results to build stochastic processes.
The proposed algorithm simulation:
1. We form geometric model of the OMU, including a platform, a output window in the fuselage and a protective radom.
2. We simulate the airflow impact on the OMU using gas dynamics for a limited time interval that is determined by the transition process and generate the output data in the form of time series of the resultant force and torque in the OMU's geometric center. The number of selected points in time series is restricted only by the required volume for storing files. Those time series are considered as some trajectories of stochastic processes. Simulation time interval is determined by the time of transients in the implementation computing circuit and computing time.
3. The obtained time series of the resultant force and torque are used to estimate spectral and correlation characteristics of stochastic processes.
4. We develop stochastic model to simulate the airflow impact on the OMU. On the basis of estimates of spectral and correlation characteristics we develop shaping filters for stochastic processes mathematical modelling.
5. We perform mathematical modelling of the airflow impact on the OMU at a predetermined time interval.
It should be noted that using of the stochastic model of the airflow on the OMU reduces the computation time by several orders compared to the gas dynamics.
Conclusion
This article describes an approach to 3D mathematical modelling of the airflow impact on the optical-mechanical unit (OMU), mounted in the aircraft's fuselage. The mathematical gas-dynamic modelling algorithm and stochastic models are used. In the proposed approach the gas dynamic model is used to retrieve the source data for constructing a stochastic model of the airflow impact on the OMU. Stochastic model is usefull for to express analysis of various design and technological solutions in the OMU's development with moderate usage of computer time. The numerical 3D gas dynamics solution is presented, airflow impact assessment on the OMU is presented.
Acknowledgements. This work is conducted with the support of grant № 16-38-60185 and of the grant № 16-01-00U4a and of the grant № 15-08-02833a.
References
1. Sobol V.R. The Mathematical Methods for Modelling of the Impact of Airflow on Optical-Mechanical Unut Based on the Results of Field Experiments. Advances of Modern Electronics, 2014, no. 4, pp. 53-57. (in Russian)
2. Goryainov A.V., Ivanov I.E., Semenikhin K.V. Identification of a Statistical Model of Atmospheric Disturbances Acting on the Optical-Dynamic Unit Placed on the Aircraft's Board During It's Flight. Abstracts of IV-th scientific and technical conference of young scientists and specialists "Actual problems of development of systems and tools ADS", Joint-stock company GSKB Almaz-Antey. 26-28 September 2013. (in Russian)
3. Weiss J.M., Maruszewski J.P., Smith W.A. Implicit Solution of Preconditioned Navier -Stokes Equations Using Algebraic Multigrid. The American Institute of Aeronautics and Astronautics Journal, 1991, vol. 37, no 1, pp. 29-36.
4. Sarkar S., Erlebacher G., Hussaini M.Y., Kreiss H.O. The Analysis and Modelling of Dilatational Terms in Compressible Turbulence. Journal of Fluid Mechanics, 1991, vol. 227, pp. 473-493.
5. Chen Y.-S., Kim S.-W. Computation of Turbulent Flows Using an Extended k-e Turbulence Closure Model. NASA Contractor Report, 179204, 1987.
6. Chieng C.C., Launder B.E. On the Calculation of Turbulent Heat Transport Downstream from an Abrupt Pipe Expansion. Numerical Heat Transfer, 1980, vol. 3, pp. 189-207.
7. Liou M.S., Steffen C.J. A New Flux Splitting Scheme. Journal of Computational Physics, 1993, vol. 107, pp. 23-39.
8. Glushko G.S., Ivanov I.E., Kryukov I.A. Method for Computation of Turbulent Supersonic Flows. Mathematical Models and Computer Simulations, 2010, vol. 2, no. 4, pp. 407-422. (in Russian)
9. Ivanov I.E., Kryukov I.A., Larina E.I. Numerical Simulation of High-Speed Viscous Intake Flow. Physico-Chemical Kinetics in Gas Dynamics, 2014, vol. 15, no. 4, 10 p. Available at: chemphys.edu.ru/issues/2014-15-4/articles/240/ (in Russian)
10. Borovikov S.N., Ivanov I.E., Kryukov I.A. Construction of Three-Dimensional Delaunay Triangulation with Constraints for Bodies with Curvilinear Geometry. Journal of Computational Mathematics and Mathematical Physics, 2005, vol. 45, no. 8, pp. 1407-1423. (in Russian)
Received February 28, 2016
УДК 517.9+519.216.2 Б01: 10.14529/ттр170411
МАТЕМАТИЧЕСКОЕ И ПРОГРАММНОЕ ОБЕСПЕЧЕНИЕ ЗБ ИМИТАЦИОННОГО МОДЕЛИРОВАНИЯ ВОЗДЕЙСТВИЯ ВОЗДУШНОГО ПОТОКА НА ОПТИКО-МЕХАНИЧЕСКОЕ УСТРОЙСТВО, РАЗМЕЩЕННОЕ В НЕГЕРМЕТИЧНОМ ОТСЕКЕ ЛЕТАТЕЛЬНОГО АППАРАТА
И.Э. Иванов, И.А. Крюков, Е.В. Ларина, В.Л. Мирошкин
Московский авиационный институт, г. Москва, Российская Федерация
Рассмотрена задача трехмерного математического моделирования воздействия воздушного потока на оптико-механический устройство (ОМУ), размещенное в негер-метизированном отсеке летательного аппарата (ЛА). Для решения указанной задачи разработана математическая модель газовой динамики, основанная на решении полной системы уравнений Навье - Стокса, которая описывают динамику турбулентного, пространственного нестационарного течения вязкого газа. Создано программное обеспечение для имитационного моделирования процесса обтекания модели ОМУ в отсеке ЛА. Воздействие воздушного потока на ОМУ описывается крутящим моментом, действующим на ОМУ со стороны воздушного потока. Приведен численный метод решения трехмерной газодинамической задачи. Численный метод базируется на численной схеме С.К. Годунова повышенного порядка точности, реализуемой на нерегулярной расчетной сетке с произвольными ячейками (тетраэдральной, призматической формы). Потоки консервативных переменных рассчитываются с помощью решения задачи Ри-мана приближенным AUSM методом. Система уравнений дополнена двухпараметриче-ской fc-моделью турбулентности, модифицированной для расчетов высокоскоростных сжимаемых течений. Для существенного сокращения затрат вычислительных ресурсов предложено использовать стохастические модели воздействия воздушного потока на ОМУ. Описывается общий алгоритм моделирования.
Ключевые слова: оптико-механическое устройство; математическое моделирование; воздействия воздушного потока.
Литература
1. Соболь, В.Р. Методика математического моделирования процесса воздействия воздушного потока на оптико-механическое устройство с учетом результатов натурных экспериментов / В.Р. Соболь // Успехи современной радиоэлектроники. - 2014. - № 4. -С. 53-57._
2. Горяинов, А.В. Идентификация статистической модели атмосферных возмущений, действующих на оптико-динамическое устройство на борту самолета во время полета / А.В. Горяинов, И.Э. Иванов, К.В. Семенихин // Тезисы докладов IV-й научно-технической конференции молодых ученых и специалистов «Актуальные вопросы развития систем и средств ВКО», ОАО ГСКБ «Алмаз-Антей». 26-28 сентября 2013.
3. Weiss, J.M. Implicit Solution of Preconditioned Navier - Stokes Equations Using Algebraic Multigrid / J.M. Weiss, J.P. Maruszewski, \Y.A. Smith // The American Institute of Aeronautics and Astronautics Journal. - 1991. - V. 37, № 1. - P. 29-36.
4. Sarkar, S. The Analysis and Modelling of Dilatational Terms in Compressible Turbulence / S. Sarkar, G. Erlebacher, M.Y. Hussaini, H.O. Kreiss // Journal of Fluid Mechanics. - 1991. V. 227. - P. 473-493.
5. Chen, Y.-S. Computation of Turbulent Flows Using an Extended k-e Turbulence Closure Model / Y.-S. Chen, S.-W. Kim // NASA Contractor Report, 179204, 1987.
6. Chieng, C.C. On the Calculation of Turbulent Heat Transport Downstream from an Abrupt Pipe Expansion / C.C. Chieng, B.E. Launder // Numerical Heat Transfer. - 1980. - V. 3. -P. 189-207. DOI: 10.1080/01495728008961754
7. Liou, M.S. A New Flux Splitting Scheme / M.S. Liou, C.J. Steffen // Journal of Computational Physics. - 1993. - V. 107. - P. 23-39.
8. Глушко, Г.С. Метод расчета турбулентных сверхзвуковых течений / Г.С. Глушко, И.Э. Иванов, И.А. Крюков // Математическое моделирование. - 2009. - Т. 21, № 12. -С. 103-121.
9. Иванов, И.Э. Численное исследование высокоскоростного течения вязкого газа в воздухозаборниках / И.Э. Иванов, И.А. Крюков, Е.В Ларина // Физико-химическая кинетика в газовой динамике. - 2014. - Т. 15, № 4. - 10 с. - URL: chemphys.edu.ru/issues/2014-15-4 / articles / 240/
10. Боровиков, С.Н. Построение трехмерной триангуляции Делоне с ограничениями для тел с криволинейной геометрией / С.Н. Боровиков, И.Э. Иванов, И.А. Крюков // Журнал вычислительной математики и математической физики. - 2005. - Т. 45, № 8. -С. 1407-1423.
Игорь Эдуардович Иванов, К&НДИДсХТ физико-математических наук, старший научный сотрудник, кафедра «Вычислительная математика и программирование:», Московский авиационный институт (г. Москва, Российская Федерация), [email protected].
Игорь Анатольевич Крюков, К&НДИДсХТ физико-математических наук, кафедра «Вычислительная математика и программирование», Московский авиационный институт (г. Москва, Российская Федерация), [email protected].
Владимировна Ларина, К&НДИДсХТ физико-математических наук, кафедра «Математическая кибернетика», Московский авиационный институт (г. Москва, Российская Федерация), [email protected].
Владимир Львович Мирошкин,
К&НДИДсХТ физико-математических наук кафедра
«» Федерация), [email protected].
Поступила в редакцию 28 февраля 2016 г.