CC BY
10www.mai.ru/science/trudy/
Труды МАИ. Выпуск № 93
УДК 629.735.45
Расчетно-экспериментальные исследования характеристик рулевого винта при вращении вертолета вокруг вертикальной оси
1 А ЛАА ЛЙЙЙ Лйййй
Леонтьев В.А. , Крымский В.С. , Игнаткин Ю.М. , Макеев П.В.
1 Центральный аэрогидродинамический институт имени Н.Е. Жуковского,
ул. Жуковского, 1, Жуковский, Московская область, 140180, Россия
Московский авиационный институт (национальный исследовательский университет), МАИ, Волоколамское шоссе, 4, Москва, A-80, ГСП-3, 125993, Россия
*e-mail: [email protected] **e-mail: [email protected] ***e-mail: [email protected] ****e-mail: [email protected]
Аннотация
В работе приведены результаты расчетно-экспериментальных исследований аэродинамических характеристик рулевого винта (РВ) при вращении вертолета вокруг вертикальной оси OYg, при отсутствии поступательной скорости. Получены и проанализированы аэродинамические характеристики изолированного РВ, а также РВ в присутствии вихревого следа несущего винта (НВ), при различных значениях угла общего шага РВ фрв и угловых скоростей вращения вертолета шу. Экспериментальные исследования выполнены на модельном ротативном стенде. Расчетные исследования проведены двумя методами: на основе дисковой квазилинейной вихревой модели винта, а также с помощью нелинейной лопастной вихревой модели винта вертолета. Полученные результаты позволяют прояснить особенности работы РВ на рассматриваемых режимах, а также способствовать
выработке рекомендаций по недопущению попадания вертолета в неуправляемое вращение.
Ключевые слова: одновинтовой вертолет, висение, вращение относительно вертикальной оси, несущий винт, рулевой винт, интерференция винтов, эксперимент, расчет.
Введение
Обеспечение безопасности полёта вертолёта является важнейшим требованием при его разработке. Опыт эксплуатации одновинтовых вертолётов с РВ при полете с малыми скоростями показывает, что при висении с боковым ветром или полете со скольжением, а также при разворотах вертолёта, РВ может попадать в область интенсивного воздействия со стороны НВ. Исследованию данного вопроса посвящен ряд экспериментальных и расчетных работ [1-5]. При развороте вертолета на режиме висения со скоростью шу, РВ может попасть в режим «вихревого кольца». Кроме того, при развороте вертолета на лопасти РВ воздействуют силы Кориолиса, вызывающие интенсивное маховое движение лопастей относительно горизонтального шарнира, что приводит к существенному изменению характера обтекания лопастей РВ.
Указанные обстоятельства могут приводить к значительному уменьшению тяги РВ, потере его эффективности, и, в результате, к неуправляемому вращению вертолёта.
Исследование данных задач требует комплексного подхода, сочетающего в
себе расчетные и экспериментальные исследования, позволяющие различными методами прояснить сущность физических процессов обтекания РВ, определить его аэродинамические характеристики, выработать рекомендации, позволяющие избежать возникновения неуправляемого вращения вертолета.
Расчетные исследования данной задачи возможны при применении современных математических моделей, учитывающих особенности формирования нелинейного вихревого следа за винтом [7; 11], а также интерференцию НВ и РВ [1; 2; 3].
В данной работе приведены результаты расчетно-экспериментальных исследований аэродинамических характеристик РВ при вращении вертолета вокруг вертикальной оси ОУё.
Экспериментальные исследования выполнены на модельном ротативном стенде. Расчетные исследования выполнены на основе математической модели движения вертолета [6], и с помощью нелинейной лопастной вихревой модели винта вертолета [7].
Экспериментальная установка
Экспериментальная установка состоит из следующих основных элементов: НВ, РВ, вращающаяся платформа, корпус вертолета, система управления и сбора информации. Общий вид экспериментальной установки без переднего обтекателя показан на рисунке 1.
V V
1 — НВ; 2 — вращающаяся платформа; 3 — корпус вертолета; 4 — РВ; Рис. 1. Общий вид экспериментальной установки
НВ имеет левое направление вращения, как это принято на отечественных вертолетах. Лопасти НВ прямоугольной формы в плане, имеют линейную крутку ДфЕ = -4°. Втулка НВ трехшарнирная, угол заклинения вала НВ вперед составляет 4°30'. Расстояние от пола лаборатории до центра втулки НВ составляет 1,8 м. Основные параметры НВ экспериментальной установки приведены в таблице 1.
Таблица 1 — Параметры НВ экспериментальной установки
Параметр Значение
Радиус НВ, Янв, мм 1250
Количество лопастей, 2л 5
Хорда лопасти, Ь, мм 87
Профиль лопасти NACA 230-10
Заполнение НВ, анв 0,107
Массовая характеристика лопасти, у 4,2
тЯнв, м/с 69
Вынос горизонтального шарнира, /гш, мм 32
Вынос вертикального шарнира, /вш, мм 73,5
Коэффициент компенсатора взмаха, к 0,42
Корпус вертолета в составе экспериментального стенда выполнен подобным фюзеляжу вертолета Ми-24 в масштабе 1:6,7.
РВ выполнен толкающим с направлением вращения «В-Н», когда верхняя лопасть при вращении винта движется назад. Расстояние от оси вала НВ до оси вала РВ /рв=1606 мм, расстояние от вертикальной плоскости симметрии модели вертолета до центра втулки РВ 2рв=150 мм. Лопасти РВ имеют прямоугольную форму в плане с профилем NACA 230-12, без крутки. Втулка РВ представляет собой совмещенный горизонтальный шарнир. Основные параметры РВ приведены в таблице 2.
Таблица 2 — Параметры РВ
Параметр Значение
Радиус РВ, Ярв, мм 293
Относительный радиус начала лопасти, г0 0,28
Количество лопастей, 2л 2
Коэффициент компенсатора взмаха, к 1
юЯрв, м/с 92
Хорда лопасти, Ь, мм 36,5
Заполнение РВ, арв 0,079
Массовая характеристика лопасти, у 0,52
Вращающаяся платформа установки служит для моделирования вращения вертолета вокруг вертикальной оси ОУё с угловой скоростью вращения шу (рисунок 1).
Программа и методика проведения эксперимента
Эксперимент проводился в приборном зале аэродинамической трубы Т-105 ЦАГИ на режимах висения с вращением, без ветра.
РВ испытывался без НВ, как изолированный, без учета индуктивного влияния вихревого следа НВ, а также в присутствии НВ с учетом индуктивного влияния вихревого следа НВ. В эксперименте без НВ значения общего шага РВ устанавливались следующими: ф7рв=1; 3,8; 6,5; 9,2; 12; 14,8; 15,9; 17; 18,2; 19,3; 20,5; 21,6°. Значения угловой скорости вращения шу изменялись в диапазоне от 0 до 132 °/с.
В эксперименте с работающим НВ значения общего шага РВ устанавливались следующими: ф7рв=14,8; 17; 19,3°. Значения угловой скорости вращения шу изменялись в диапазоне от 0 до 90 °/с. Общий шаг НВ для всех режимов устанавливался на величину ф7нв=7° Коэффициент тяги НВ при этом составляет Ст =0,0107, а коэффициент крутящего момента шк=0,0009.
Во время эксперимента при помощи регистрирующей аппаратуры записывались величины тяги РВ ТРВ, крутящего момента РВ МкРВ, а также угловой скорости вращения вертолета шу.
Описание математической модели движения вертолета FDH_tr ЦАГИ
Математическая модель FDH_tr позволяет проводить расчет неустановившегося движения вертолета. При этом аэродинамические характеристики РВ и НВ определяются интегрированием погонных сил и момента, действующих в сечениях лопастей по радиусу и азимутальному положению лопасти. Подробное описание уравнений движения упругих лопастей и их решение содержится в работе [6]. Аэродинамические нагрузки на лопасть определяются на основе гипотезы плоских сечений, с учетом изменения угла атаки а и чисел Маха М в диапазоне 0<а<360°, 0<М<1. Коэффициенты аэродинамических характеристик профилей, для заданных чисел М и Яв определяются путем обработки нелинейных экспериментальных зависимостей, полученных при испытаниях профилей в аэродинамических трубах. Для вычисления скоростей, индуцируемых винтом в собственной плоскости вращения, используется квазилинейная дисковая вихревая теория винта Е.С. Вождаева [8], в диапазоне скоростей обдувки винта, когда вихревой след можно считать линейным.
Описание нелинейной лопастной вихревой модели винта
В основе нелинейной лопастной вихревой модели винта [7; 9] лежит теория несущей линии. Для расчета аэродинамических характеристик лопасти используется гипотеза плоских сечений. Аэродинамические характеристики профилей определялись с использованием данных экспериментальных продувок при соответствующих значениях чисел Рейнольдса и Маха. В расчетах учитывается маховое движение лопастей. Свободные вихри, отходящие от лопастей,
располагаются по линиям тока. В ходе расчета за винтом выстраивается пространственный нелинейный вихревой след. Особенностью используемой в модели вихревой системы, является представление вихревых отрезков пелены в виде диффундирующих вихревых линий, что позволяет более точно моделировать физические процессы, происходящие в вихревой пелене.
В данных расчетах была использована расчетная модель комбинации НВ и РВ, с учетом их интерференции, изображенная на рисунке 2.
Вихревой след за лопастью НВ
Рис. 2. Расчетная схема комбинации НВ и РВ стенда Экспериментальные и расчетные поляры РВ на режиме висения, без влияния
НВ
На рисунке 3 показаны экспериментальные и расчетные поляры РВ без влияния индуктивного потока от НВ. Видно, что до угла установки ф7рв ~ 15° поляры РВ, полученные на основе расчетов удовлетворительно согласуются с
экспериментальной полярой. Далее наблюдаются отличия в полярах, но начиная с ф7рв ~ 18° все три поляры вновь близки.
0,001 0,002 0,003 0,004 0,005
Рис. 3. Экспериментальные и расчетные поляры РВ на режиме висения, без влияния
НВ
Влияние сил Кориолиса на маховое движение лопастей РВ при вращении вертолета относительно вертикальной оси
Особенностью работы РВ при вращении вертолета относительно оси OYg, с угловой скоростью шу, является возникновение на лопастях РВ сил Кориолиса Fk. Это происходит в силу того, что лопасти РВ, вращающиеся с угловой скоростью вокруг оси РВ - юРВу вращаются также с угловой скоростью шу (см. рис. 4).
Трв
Рис. 4. Схема возникновения и направление действия сил Кориолиса Бк на лопастях шарнирного РВ при вращении вертолета относительно вертикальной оси со
скоростью юу
На рисунке 4 условно показано направления действия сил Кориолиса при направлении вращения РВ «В-Н» (верхняя лопасть движется назад). Видно, что направление сил Кориолиса различно для нижней лопасти, движущейся в сторону оси вращения вертолета, и верхней лопасти, движущейся от оси. Величина сил Кориолиса зависит от произведения угловых скоростей юрв-юу и с увеличением шу возрастает, вызывая более интенсивное маховое движение лопастей РВ.
На рисунке 5 для одной из лопастей РВ показано изменение угла взмаха лопасти вл по углу азимута щ. Видно, что в данном случае, при определенной величине юу, максимальный и минимальный углы взмаха лопасти достигают предельно возможных значений, устанавливаемых конструктивными ограничителями взмаха лопастей.
90
Рис. 5. Маховое движение лопасти при вращении вертолета относительно вертикальной оси ОУё, со скоростью юу
При наличии интенсивного махового движения лопастей возможен срыв потока на лопасти, отклоняющейся вниз от плоскости вращения РВ. На другой части диска РВ, где лопасть отклоняется, соответственно, вверх, наблюдается уменьшение угла атаки лопастей, что в результате может привести к падению тяги РВ [10]. Также большое значение при этом играет величина компенсатора взмаха лопастей РВ, снижающего угол установки лопасти при ее взмахе вверх и, напротив увеличивающего угол установки при взмахе лопасти вниз, что может усилить эффект уменьшения тяги РВ при угловом вращении вертолета со скоростью шу.
Экспериментальные и расчетные аэродинамические характеристики РВ при вращении вертолета относительно вертикальной оси, без НВ
На рисунках 6-13 представлены экспериментальные и расчетные зависимости коэффициентов тяги и крутящего момента РВ Ст =/(юу) и ти=/(юу) для ф7рв=3,8; 9,2; 14,8; 21,6°.
Рис. 6. Экспериментальная и расчетные зависимости Ст =/(юу) РВ, при ф7рв=3,8°
тк
0,002 0,0015 0,001
0,0005 $ 0
II II «Р7рв= =3,8°
Л рт
-й—Расчет МАИ (изолированный РВ)
— - Л
£г ■п.
- — — —
0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 СОу,%
Рис. 7. Экспериментальная и расчетная зависимости шк=/(юу) РВ, при ф7рв=3,8°
Рис. 8. Экспериментальная и расчетные зависимости Ст =/(<лу) РВ, при ф7рв=9,2°
тк
0,0025 0,002 0,0015 0,001 0,0005 0
V- £— - -- Г (Р1р 2°
- -К - 3- В
о— Эксперимент ЦАГИ (изолированный РВ) ¿г- Расчет МАИ (изолированный РВ) —
-1-1-1-
10 20 30 40 50 60 70 80 90 100 110
120 130 СО у, %
Рис. 9. Экспериментальная и расчетная зависимости тк=/(юу) РВ, при ф7рв=9,2°
Рис. 10. Экспериментальная и расчетные зависимости Ст =/(юу) РВ, при ф7рв=14,8°
тк
0,0035 0,003 0,0025 0,002 0,0015 0,001 0,0005 О
■л
-с у
--
——С -— >в 5
<р 14
Л 1 1 ЛГ1Л / - ПО\
-а—Расчет МАИ (изолиоованный РВ)
—1 — —1 г-
10 20 30 40 50 60 70 80 90 100 110
120 130 СОу, %
Рис. 11. Экспериментальная и расчетная зависимости тк=/(шу) РВ, при ф7рв=14,8°
Рис. 12 Экспериментальная и расчетные зависимости Ст =/(юу) РВ, при ф7рв=19,3°
Рис. 13. Экспериментальная и расчетная зависимости шк=/(шу) РВ, при ф7рв=19,3°
Из представленных зависимостей можно заключить, что в целом наблюдается удовлетворительное соответствие результатов расчетов с экспериментальными данными. Имеющиеся отличия расчетных и экспериментальных зависимостей можно объяснить принятыми в математических моделях допущениями.
Из графиков зависимостей Ст =/(ыу) следует, что при малых углах установки лопастей РВ ф7рв=3,8° (рис. 6), с ростом угловой скорости вращения юу, тяга РВ возрастает, а начиная с ф7рв=9,2° и до максимального, с ростом угловой скорости вращения вертолета шу тяга РВ начинает снижаться. Результаты экспериментальных исследований подтверждают результаты анализа и математического моделирования,
проведенных в работе [10]. Снижение тяги РВ связано не столько с режимами
14
«вихревого кольца» РВ, хотя они также имеют место при больших величинах угловых скоростей, сколько с интенсивным маховым движением лопастей, возникающим из-за сил Кориолиса, как это было отмечено выше.
На рисунке 14 представлены диаграммы распределения истинного угла атаки элемента лопасти а по диску РВ, полученные на основе обработки расчетов.
Диаграммы построены для ф7рв=14,8° и юу=0; 25; 50; 75; 100; 125 °/с с помощью цветового градиента, отражающего интенсивность рассматриваемого параметра, и представлены для одного оборота РВ.
Из приведенных на рисунке 14 диаграмм распределения по диску РВ угла атаки сечения лопасти видно, что на верхней части диска с ростом скорости вращения шу и при интенсивном маховом движении лопастей, наблюдаются минимальные углы атаки в сечениях лопасти, близкие к нулю, а на нижней части диска, напротив, чрезвычайно большие величины угла атаки, лежащие за срывными зонами. Данное обстоятельство связано с интенсивным маховым движением лопастей с большой амплитудой угла взмаха лопасти, а также с наличием большой величины компенсатора взмаха лопасти к=1 для рассмотренной втулки РВ ротативного стенда с совмещенным ГШ. На верхней части диска РВ лопасть взмахивает вверх и за счет компенсатора взмаха происходит уменьшение угла ее установки, а на нижней части лопасть движется вниз, что приводит к увеличению угла установки и достижению критических углов атаки в сечениях лопасти. Указанные причины приводят к резкому падению тяги практически до нуля на
верхней части диска РВ и значительному снижению тяги РВ на нижней части диска РВ из-за срыва потока на лопасти.
(Оу=0 % (Оу=25 7с
Рис. 14. Диаграммы распределения углов атаки а сечения лопасти по диску РВ для
ф7рв=14,8° и Юу=0; 25; 50; 75; 100; 125 °/с
Экспериментальные и расчетные аэродинамические характеристики РВ при вращении вертолета относительно вертикальной оси совместно с НВ
На рисунках 15, 16 представлены экспериментальные и расчетные зависимости коэффициентов тяги и крутящего момента РВ Ст =Дюу) и тк=/(юу), полученные с учетом индуктивного влияния вихревого следа НВ, для ф7рв=14,8° и ф7нв=7°. Для сравнения там же приведены зависимости для изолированного РВ, которые были рассмотрены выше.
Из приведенных зависимостей можно сделать вывод о том, что индуктивное влияние вихревого следа НВ на аэродинамические характеристики РВ незначительно (для случая висения, когда вертолет вращается относительно вертикальной оси на месте без поступательной скорости).
-
г/ >7п = 1 4, а 7 7 О — ■
■Расчет ЦАГИ (изолированный Эксперимент ЦАГИ (изопиров р Я) ЙР ▲ -
аннь 3)
■Расчет МАИ (изолированный РВ) Эксперимент ЦАГИ (с учетом рабо -Расчет МАИ (с учетом работы НВ)
ты НВ)
—I — — — 1 —I- —г —г —г
О 10 20 30 40 50 60 70 80 90 100 110 120 130 СОг,%
Рис. 15. Зависимости Ст =Дюу) изолированного РВ и РВ в присутствии вихревого
следа НВ, при ф7рв=14,8°
тк
0,0035 0,003 0,0025 0,002 0,0015 0,001 0,0005 0
— г Д
-- — -- ■А
7 4 ЯЯ гг V
А— исг 1
(р7 = 1 4,8 о. Ч>1 --Т
рв НВ
—□—Эксперимент ЦАГИ (изолированный РВ) а Расчет МАИ (изолированный РВ) —
—•—Э (сперимент ЦАГИ (с учетом работы НВ) эсчет МАИ (с учетом работы НВ) —
—тк—Рг
10 20 30 40 50 60 70 80 90 100 110 120
соу, °/с
Рис. 16. Зависимости тк=/(шу) изолированного РВ и РВ в присутствии вихревого
следа НВ, при ф7рв=14,8°
На рисунке 17 приведена расчетная форма вихревого следа НВ и РВ для одной из расчетных точек, приведенной на рисунках 15 и 16 - для шу=100 °/с. Можно отметить сложную и существенно нелинейную картину вихревого следа РВ, свидетельствующую о неравномерности его обтекания, которая происходит из-за интенсивного махового движения лопастей, вызванного силами Кориолиса. При этом, верхняя часть диска ванта, где углы атаки сечений лопасти минимальны, работает на режимах близких в авторотации, а нижняя - где углы велики и имеют срывные значения - на режиме близком к «вихревому кольцу».
Рис. 17. Визуализация расчетной формы вихревого следа РВ с учетом индуктивного воздействия вихревого следа НВ, выполненная по модели МАИ, для ф7рв =14,8° и
а>=100 °/с
Заключение
Проведено сравнение результатов расчета аэродинамических характеристик РВ модельного ротативного стенда на базе двух математических моделей.
Сравнение расчетной поляры изолированного РВ с данными эксперимента, показало удовлетворительное соответствие.
Приведены результаты расчета аэродинамических характеристик изолированного РВ при вращении с угловыми скоростями шу=0...125 °/с, для фиксированных углов установки лопастей ф7рв=3,8; 9,2; 14,8; 19,3°.
Сравнение экспериментальных зависимостей Ст=/(шу) и Ш1=Дюу) с расчетными на режимах вращения показывает, в целом, удовлетворительное соответствие. Имеющиеся отличия расчетных и экспериментальных зависимостей можно объяснить принятыми в математических моделях допущениями.
В результате анализа полученных данных установлено, что с увеличением окружной скорости вращения шу на лопасти РВ, кроме аэродинамических сил, действуют значительные нагрузки, связанные с возникновением силы Кориолиса, имеющие разное направление на верхней и нижней частях диска РВ, вызывающие интенсивное маховое движение лопастей РВ. Интенсивность махового движения увеличивается с ростом скорости вращения шу.
Из приведенных диаграмм распределения по диску РВ угла атаки сечения лопасти видно, что на верхней части диска с ростом скорости вращения шу наблюдаются минимальные углы атаки в сечениях лопасти, близкие к нулю, а на нижней части диска, напротив, чрезвычайно большие величины угла атаки, лежащие в зоне срыва. Указанные причины приводят на средних и больших ф7рв к резкому падению тяги практически до нуля на верхней части диска РВ и значительному снижению тяги РВ на нижней части диска РВ.
Результаты экспериментальных и расчетных исследований показывают, что при небольших ф7рв с увеличением шу, происходит существенное увеличение тяги
РВ. Сечения лопасти РВ находятся на досрывных углах атаки. Для средних значений ф7рв эффект от роста скорости вращения вертолета шу проявляется иначе. Так при ф7рв=9,2° с ростом скорости вращения шу тяга РВ меняется слабо, но при достижении скорости шу=100 °/с, начинает снижаться. Сечения лопасти РВ начинают попадать в зону срыва потока. При больших ф7рв=14,8°; 19,3°, с ростом скорости вращения шу, наблюдается существенное снижения тяги РВ. Сечения лопасти РВ находятся в зоне срыва.
Расчетные зависимости Сг=/(шу) и тк=/(шу) РВ на режимах вращения вертолета совместно с НВ, удовлетворительно согласуются с экспериментом и показывают, что вихревой след НВ незначительно влияет на аэродинамические характеристики РВ.
Таким образом, сравнение расчетных исследований с экспериментальными данными показывает, что разработанные расчетные модели могут быть использованы при исследовании аэродинамических характеристик РВ на режиме вращения вертолета относительно вертикальной оси.
Библиографический список
1. Игнаткин Ю.М., Макеев П.В, Шомов А.И. Интерференция несущего и рулевого винтов вертолета при полете со скольжением // Труды МАИ. 2015. №82 ШЬ: http://www.mai.ru/science/trudy/published.php?ID=58605
2. Игнаткин Ю.М., Макеев П.В., Шомов А.И. Численное моделирование
интерференции несущего и рулевого винтов вертолета на базе нелинейной
21
лопастной вихревой модели // Известия высших учебных заведений. Авиационная техника. 2016. № 1. С. 41-47.
3. Щеглова BM. Bихpевая структура комбинации несущий винт - рулевой винт на большой скорости полета // Научный вестник МГТУ ГА. 2013. №188. С 137142.
4. Бутов B.n. Структура, геометрия и интенсивность спутного вихревого следа несущих винтов одновинтовых и соосных вертолетов в реальных условиях полета // Труды 4-го форума Pоссийского вертолетного общества. Москва, 2000. C 19 - 34.
5. Timothy M. Fletcher and Richard Brown. Helicopter tail rotor thrust and main rotor wake coupling in crosswind flight // Journal of Aircraft. 2010. №47 (6). pp. 21362148.
6. Леонтьев BA. Метод решения уравнений движения упругих лопастей вертолетных винтов в общем случае движения // Ученые записки ЦАГИ. Т. XLI. № 5. 2010. С 67-79.
7. Игнаткин Ю.М., Макеев n.B., Гревцов Б.С., Шомов А.И. Нелинейная лопастная вихревая теория винта и ее приложения для расчета аэродинамических характеристик несущих и рулевых винтов вертолета // Bестник Московского авиационного института. 2009. Т.16. №5. С. 24-31.0
8. Bождаев Е.С. К теории несущего винта с произвольным углом атаки // Ученые записки ЦАГИ. 1972. Т.3 №2. С. 1-10.
9. Игнаткин Ю.М., Макеев n.B., Шомов А.И. Численное моделирование прикладных задач аэродинамики вертолета на базе нелинейной лопастной
вихревой модели винта // Труды МАИ. 2016. №87. URL:
http://www.mai.ru/science/trudy/published.php?ID=58605
10. Анимица В.А., Леонтьев В. А. О «самопроизвольном» вращении одновинтовых вертолетов // Научный вестник МГТУ ГА. 2010. №151. С. 96-102.
11. Вершков В.А., Крицкий Б.С., Махнев М.С., Миргазов Р.М., Требунских Т.В. Сравнение результатов численного моделирования обтекания несущего винта в различных пакетах программ // Труды МАИ, 2016. №89. URL: http://www.mai.ru/science/trudy/published.php?ID=72704
