CC BY
13УДК 533.6.011.8:004.94
моделирование сил и моментов Сил набегающего потока атмосферы в целях верификации динамических режимов системы управления движением и навигации мкс и синтеза оптимального управления
© 2017 г. Атрошенков С.н.1, прутько A.A.1, Крылов А.н.1, Крылов н.А.1, Тубарев Ф.в.2
'Ракетно-космическая корпорация «Энергия» им. С.П. Королёва» (РКК «Энергия») Ул. Ленина, 4А, г. Королёв, Московская обл., Российская Федерация, 141070, e-mail: [email protected]
2ООО «ДАТАДВАНС» Научный проезд, 17, г. Москва, Российская Федерация, 117246, e-mail: [email protected]
Рассматриваются два метода моделирования сил и моментов сил от воздействия набегающего потока атмосферы на поверхность Международной космической станции (МКС) в процессе выполнения динамических режимов системы управления движением и навигации служебного модуля, которые разработаны РКК «Энергия». Эти методы, пригодные как для поиска оптимального управления движением МКС, так и для верификации динамических режимов системы управления движением и навигации служебного модуля, основаны на расчете аэродинамических характеристик в наборе точек четырехмерного (4D) пространства и последующей 4D-uнmepnoляцuu аэродинамических характеристик в ходе моделирования. Сравниваются результаты расчетов воздействия набегающего потока атмосферы во время моделирования оптимального по расходу топлива разворота (optimal propellant maneuvers) МКС на большой угол. Во время разворота геометрия МКС существенно меняется за счет поворотов крупногабаритных элементов МКС вокруг разнонаправленных осей. Предложенные методы планируется применить при подготовке космического эксперимента «МКС-Разворот». Они также могут быть пригодны при разработке алгоритмов управления системы управления движением и навигации будущих крупногабаритных космических объектов.
Ключевые слова: аэродинамические характеристики, четырехмерное (4D) пространство, 4D-интерполяция, оптимальный разворот (optimal propellant maneuvers), методы вычисления аэродинамических характеристик.
AERODYNAMIC FORCES AND TORQUES SIMULATION FOR THE VERIFICATION OF INTERNATIONAL SPACE STATION GUIDANCE, NAVIGATION AND CONTROL SYSTEM DYNAMIC MODES
AND FOR OpTIMAL CONTROL SYNTHESIS Atroshenkov S.N.1, Prutko A.A.1, Krylov A.N.1, Krylov N.A.1, Gubarev F.V.2
1S.P. Korolev Rocket and Space Public Corporation Energia (RSC Energia) 4A Lenin str., Korolev, Moscow region, 141070, Russian Federation, e-mail:[email protected]
2DATADVANCE LLC
17 Nauchny proezd, Moscow, 117246, Russian Federation, e-mail: [email protected]
The paper reviews two methods of simulating forces and torques produced by the flow of atmospheric air incident on the surface of the International Space Station (ISS) during execution of dynamic modes of the Guidance, Navigation and Control System of the Service Module, which were developed by RSC Energia. These methods suitable both for searching for the optimal ISS motion control, and for the verification of dynamic modes of the Service Module Guidance, Navigation and Control System, are based on the computation of aerodynamic properties within set of points in a 4D-space and subsequent 4D interpolation of aerodynamic properties during simulation. The paper compares the results of computations of exposure to incident air flow during simulations of the ISS optimal-propellant wide-angle turning maneuver.
During a turning maneuver the ISS geometry changes significantly due to large-scale elements of the ISS rotating about axes pointing in various directions. There are plans to use the proposed methods during preparations for space experiment MKS-Razvorot. They can also be useful in developing control algorithms for the Guidance, Navigation and Control Systems of future large-scale space vehicles.
Key words: aerodynamic characteristics, 4D-space, 4D-interpolation, optimal propellant maneuver, methods for computing aerodynamic characteristics.
АТРОШЕНКОВ С.Н.
ПРУТЬКО А.А.
КРЫЛОВ А.Н.
КРЫЛОВ Н.А.
ГУБАРЕВ Ф.В.
АТРОШЕНКОВ Сергей Николаевич — ведущий инженер-математик РКК Энергия, e-mail: [email protected]
ATROSHENKOV Sergey Nikolaevich — Lead engineer-mathematician at RSC Energia, e-mail: [email protected]
ПРУТЬКО Алексей Александрович — инженер-математик 2 категории РКК «Энергия», e-mail: [email protected]
PRUTKO Aleksey Aleksandrovich — Engineer-mathematician 2 category at RSC Energia, e-mail: [email protected]
КРЫЛОВ Андрей Николаевич — кандидат физико-математических наук, начальник сектора РКК «Энергия», e-mail: [email protected]
KRYLOV Andrey Nikolaevich — Candidate of Science (Physics and Mathematics), Head of Subdepartment at RSC Energia, e-mail: [email protected]
КРЫЛОВ Николай Андреевич — инженер 2-й категории РКК «Энергия», e-mail: [email protected]
KRYLOV Nikolay Andreevich — Engineer 2 category at RSC Energia, e-mail: [email protected]
ГУБАРЕВ Федор Васильевич — начальник отдела OOO «ДАТАДВАНС», e-mail: [email protected]
GUBAREV Fedor Vasilyevich — Head of Department at DATADVANCE LLC, e-mail: [email protected]
введение
В процессе полета Международной космической станции (МКС) бортовая система управления движением и навигации (СУДН) выполняет различные динамические режимы: программные развороты, поддержание заданной ориентации, стабилизация вектора тяги во время коррекции орбиты, неуправляемый полет в индикаторном режиме и т. д. Верификация каждого такого режима и их совокупности на длительных интервалах полета МКС требует проведения предварительного моделирования как для верификации работоспособности СУДН, так и для планирования потребных ресурсов, в частности — расхода топлива для ракетных двигателей (РД) объединенной двигательной установки МКС. Необходимо моделировать основные внешние факторы, действующие на МКС — силы и моменты сил, создаваемые гравитацией («гравитационные»), а также силы и моменты, создаваемые набегающим потоком атмосферы (далее — «аэродинамические сила и момент»). В настоящее время МКС — объект с непрерывно меняющейся геометрией за счет вращения по нескольким степеням свободы панелей солнечных батарей (СБ) и радиаторов Американского сегмента (АС) МКС. При этом положение Солнца отслеживают не все поворотные элементы. Так, на участке полета МКС при выполнении операции сближения с российским космическим аппаратом американские коллеги фиксируют приводы «Альфа» своих СБ и позволяют отслеживать Солнце приводами СБ «Бета», при этом некоторые приводы «Бета» фиксируют. Схематически расположение поворотных элементов МКС иллюстрирует рис. 1. Описание элементов МКС содержится в работе [1].
Учитывая сказанное, понятно, что расчет аэродинамических силовых факторов, действующих на МКС, геометрия которой непрерывно меняется во время выполнения динамического режима, является непростой задачей.
Рис. 1. Вид МКС с поворотными элементами: 1В, 3В, 1А,
3А — приводы «Бета» солнечных батарей (СБ) Американского сегмента (АС) правого борта; 2А, 4А, 2В, 4В — приводы СБ АС левого борта; а_1 и а_2 — приводы «Альфа» правого и левого бортов, соответственно; Т_1 и Т_2 — приводы большихрадиато-ров правого и левого бортов, соответственно Примечание. Описание элементов МКС содержится в работе [1]. Здесь ориентация МКС относительно осей орбитальной системы координат произвольна.
Для моделирования текущего аэродинамического момента, действующего на МКС, РКК «Энергия» применяла различные методики в зависимости от стадии развертывания МКС. В данной работе описан метод, разработанный для последних стадий развертывания МКС.
1. метод расчета аэродинамических характеристик в ходе моделирования
Выражения для векторов результирующей силы Fа и главного момента Ta аэродинамического сопротивления имеют вид [2, 3]:
F = Q CF (0; T = Q ^,(0, (1)
где ^^ = 0,5р(О(70(О)2 — скоростной напор набегающего потока. Скалярные величины р и У0 — это плотность набегающего потока газа и его скорость относительно МКС. Векторы CF и CT — векторы аэродинамических характеристик (АДХ), которые зависят от геометрии МКС
и направления набегающего потока. Параметр V0(t) вычисляется в модели движения центра масс (ЦМ) МКС по околоземной орбите. Плотность атмосферы p(t) вычисляется с помощью моделей движения ЦМ МКС и Солнца, а также в соответствии с ГОСТ [4] и по прогнозам солнечной и геомагнитной активности Центра Маршалла NASA. Расчет скалярных параметров не представляет особых вычислительных трудностей.
В то же время, расчет АДХ (сил CF и моментов CT) — это трудоемкая вычислительная задача даже для одного, фиксированного момента времени. В соответствии с методиками, принятыми РКК «Энергия», поверхность МКС текущей геометрии разбивается на большое количество N плоских элементов, и вычисляются аэродинамические силы CF действующие на каждый элемент n, и точка ее приложения гм внутри элемента. При вычислении CF учитываются молекулярный состав набегающего потока, свойства материала элемента и его температура, поток молекул, достигающих его поверхности из-за эффектов затенения другими элементами. Вычисления проводятся в системе координат (СК), связанной с корпусом МКС для одного (текущего) направления набегающего потока a . После вычисления всех
^ aer
аэродинамических коэффициентов сил для всех N элементов текущие АДХ МКС задаются выражениями:
N N N
CF = X CF; Ct = X C1= X [г,- x CF]. (2)
i - 1 i-1 i -1
Такой способ вычисления АДХ далее в тексте именуется «прямой».
Число N панелей разбиения для современной конфигурации МКС составляет несколько сотен тысяч, и вычисление АДХ требует больших вычислительных ресурсов даже на современных персональных ЭВМ и кластерах. Так, расчет АДХ для одного варианта геометрии МКС и одного направления потока требует нескольких секунд. Полученные АДХ имеют погрешность: ACf = 15% — для сил и АСТ = 30% — для моментов. Моделирование движения МКС в ходе разработки и верификации бортовых алгоритмов управления требует вычислений в режиме реального времени либо еще быстрее, и в таких условиях реализовать прямой расчет АДХ не всегда возможно. Для вычисления АДХ в ходе моделирования движения объекта с изменяющейся геометрией требуются подходы, отличные от описанных выше.
Для расчета аэродинамических сил и моментов в ходе моделирования динамических режимов МКС РКК «Энергия» разработан многоступенчатый процесс, в котором участвуют специалисты по бортовому программному обеспечению СУДН, специалисты-расчетчики аэродинамических характеристик, а также специалисты ЦУП-М управления полетом МКС.
Идея метода основана на предположении, что в каждый момент времени выполнения моделируемого динамического режима МКС нам известен закон, по которому вращаются все панели СБ и радиаторы МКС. В частности, известно, в каких положениях зафиксированы приводы «Бета» неподвижных СБ АС и какой режим вращения задан для остальных приводов «Бета» СБ АС: отслеживание Солнца либо отслеживание Солнца с постоянным угловым смещением. Как правило, такая информация выдается американскими коллегами за месяц до исследуемой динамической операции. В таком случае текущие положения всех поворотных элементов и текущая геометрия МКС могут быть рассчитаны в зависимости от единичного вектора еу текущего направления на Солнце в СК, связанной с корпусом МКС. В базовой СК служебного модуля (СМ) е5 задается двумя углами: углом vС местного склонения Солнца, т. е. углом между плоскостью ХСМ1СМ и направлением на Солнце, и углом 9С местного азимута Солнца, который отсчитывается в плоскости ХСМ^СМ между осью ХСМ и проекцией направления на Солнце на плоскость ХСМ2СМ. Единичный вектор ааег также определяется двумя углами: углом атаки аа и углом скольжения Ра [2]. Таким образом, четверка углов {аа, ва, vС, 6С| полностью определяет исходные данные для расчета АДХ (рис. 2).
Рис. 2. Углы-параметры 40-пространства в системе координат служебного модуля МКС. Здесь угол атаки аа отрицателен
Далее проводится предварительный анализ движения МКС в исследуемом динамическом режиме, и определяются пределы изменений каждого из четырех углов. Затем в четырехмерном (4Б) пространстве проводится конечно-элементное (КЭ) разбиение областей изменения Е(аа, Ра, vС, 0С) угловых параметров на 4Б-«прямоугольные» призмы так, чтобы объединение Е(аа, ра, ус, 9с) КЭ-призм содержало Е(аа, ра, ус, 0С). Обозначим К4(аа, Ра, ус, 0С) множество узлов-вершин в КЭ-разбиении Е(аа, ра, ус, 0с). Процесс вычисления АДХ для МКС с меняющейся геометрией выглядит следующим образом.
На первом, предварительном этапе вычисляются векторы {с , с , с , m , m , m }
1 Уху г X у г>
(в дальнейшем — «6-векторы») АДХ для каждого узла из К4(аа, Ра, ус, 0С). Это делается с помощью следующей процедуры. Пусть область изменения для направлений набегающего потока {аа, Ра} в К4 состоит из набора 2Б-узлов {а'а, р а }, I = 1, 2, ..., I; ] = 1, 2, ..., ], а область изменения параметров положения Солнца {у^ 0с} в К4 состоит из набора 2Б-узлов {уС, 0С }, где к = 1, 2, ..., К; I = 1, 2, ..., Ь. Для каждого сочетания (к, I), зная законы вращения/ фиксации поворотных элементов МКС в зависимости от положения Солнца, вычисляем набор Бк1{йр} положений йр всех приводов поворотных элементов МКС. На данном этапе развертывания МКС функционируют 16 приводов вращения, т. е. р = 1, 2, ..., 16. Обозначим итоговый результат расчета как К2{ус, 0С, Б} — множество из КЬ строк по 16 элементов в строке. По этим данным рассчитываются «узловые» АДХ для каждого из узлов КЭ-разбиения Е(аа, ра, ус, 0с), причем в процессе расчета для каждой строки Бы из К2{ус, 0С, Б} производятся следующие действия:
• построение геометрии МКС и ее разбиение на элементарные панели;
• расчет I-] 6-векторов АДХ для каждого 2Б-узла {а!а, р а} из области изменения {а , Р } множества К4.
аа
В результате в каждом узле КЭ-разбие-ния Е(аа, Ра, ус, 0С) вычисляется 6-вектор (С^, Cтт)т. Совокупность АДХ для всех узлов Е(аа, Ра, vС, 0С) назовем 4Б-базой В4Б.
Второй этап — это вычисление АДХ в ходе моделирования динамического режима. Схема вычислений для текущего момента времени моделирования t такова:
• вычислить «положение» в 4Б-пространстве Ъ(0 = {аа(0, Ра(0, vС(í), 0С(О};
• найти КЭ-призму Мст из разбиения Е такую, что точка Ъ(0 лежит либо внутри призмы Мсик, либо на ее границе;
• найти в 4Б-базе 16 наборов АДХ (узловые АДХ) для 16-ти узлов-вершин 4Б-призмы Ыси
• вычислить текущие АДХ по значениям узловых АДХ и текущим значениям Ъ(0 с помощью интерполяции. В данной работе представлены результаты интерполяции АДХ внутри четырехмерного конечного элемента первого порядка с помощью Сирендипова семейства многочленов путем обобщения на четырехмерное пространство техники, описанной, например, в работе [5].
Предварительный анализ и минимизация числа узлов в КЭ-разбиении Е(аа, Ра, vС, 0С) является отдельной и очень важной задачей. В качестве примера рассмотрим КЭ-разбиение Е10, покрывающее всю область определения 4Б-пространства {аа, Ра, ус, 0С}. Пусть каждая призма — 4Б-куб со стороной 10°. Тогда по оси аа имеем I = 37 узлов (-180 < аа <180°); по оси Ра — ] = 19 узлов; по оси vСfl— К = 19 узлов (-90° < Ра; vС < 90°); по оси 0С — Ь = 37 узлов (0 < 0С < 360°). Число 4Б-кубов в КЭ-разбиении Е10 равно 419 904, общее число узлов 494 209. Если принять, что расчет АДХ стандартного качества для одного узла требует 3 с, то для расчета 4Б-базы потребуется более 400 ч процессорного времени. Очевидно, что для расчета 4Б-базы желательно выделить минимальное подмножество Ет.п из Е10, содержащее траекторию движения МКС для моделируемого динамического режима в 4Б-пространстве.
Для некоторых динамических режимов, таких как стабилизация относительно орбитальной системы координат (ОСК) или инерциальной системы координат (ИСК), разворот на небольшой угол, непродолжительный неуправляемый полет и др., выделить подмножество Ет.п сравнительно несложно.
Более сложный случай возникает при синтезе оптимальных разворотов МКС на большие углы относительно ОСК, когда оптимальную траекторию разворота предсказать очень трудно.
2. оптимальный разворот мкС на 180° вокруг оси УоСк
В качестве иллюстрации рассмотрим верификацию разворота МКС относительно оси УОСК на угол =180° длительностью
Т = 5 500 с. Оси ОСК определены слешах
дующим образом: ось YOCK направлена вдоль радиуса-вектора из центра Земли в ЦМ МКС, ось ZOCK направлена против вектора орбитальной угловой скорости движения ЦМ МКС, ось ХОСК дополняет тройку до правой и в случае круговой орбиты направлена вдоль вектора линейной скорости ЦМ МКС.
Ориентацию базовой системы координат СМ XYZCM относительно ОСК задаем с помощью собственного кватерниона Q6. Если оси XYZCM базовой СК СМ совпадают с осями ОСК, то говорят, что имеет место ориентация МКС «ОСК» (в этом случае qb = (1, 0, 0, 0)). Если оси базовой СК СМ совпадают с осями ОСК, развернутой на 180° вокруг оси Y^^ то говорят, что имеет место ориентация МКС «ОСК-1» (в этом случае Q6 = (0, 0, 1, 0)). Во время полета МКС подолгу поддерживает ориентацию в положениях TEA (Torque Equilibrium Attitude), в которых сумма действующих на МКС моментов сил гравитации, аэродинамики и инерции равна нулю. Рассматриваемый разворот начинается из положения ТЕА, которое близко к ориентации ОСК-1, и заканчивается в положении ОСК. При этом угловые скорости МКС относительно ОСК в начальный и конечный моменты времени равны нулю. Такой разворот NASA именует «разворот mXVV-pXVV». Развороты такого типа, а также обратные развороты, в настоящее время требуются перед выполнением операции стыковки или расстыковки российских транспортных модулей.
Траектории оптимальных по расходу топлива разворотов МКС с такими начальными и конечными условиями были первоначально рассчитаны подрядчиком NASA — фирмой Charles Stark Draper Laboratories и именуются в литературе Optimal Propellant Maneuvers (OPM) [6-8]. Далее в тексте эти развороты и их траектории именуются ОРМ NASA. Заметим, что задача оптимального управления по расчету ОРМ решалась для объекта управления — твердого тела; также компоненты вектора управления трактовались как гладкие функции, лишь по максимальной величине соответствующие моментам сил, создаваемым РД. Длительность разворота фиксирована и составляет 5 500 с. Заметим также, что задача расчета ОРМ NASA была решена численно, при помощи коммерческого программного продукта DIDO с интерфейсом системы MATLAB [6, 7].
Начиная с 2012 г., МКС регулярно выполняет развороты под управлением АС МКС, используя при этом РД Российского сегмента (РС) МКС в режиме USTO (US Thrusters only). Режим USTO был включен в состав СУДН МКС по просьбе NASA, чтобы иметь возможность выполнять разгрузку кинетического момента силовых гиродинов (CMG) под управлением АС МКС. Разворот МКС под управлением АС выполняется путем отслеживания оптимальной траектории ОРМ NASA с помощью включений РД РС по циклограмме разгрузки CMG.
Российские специалисты также проводят работы по оптимальному по расходу топлива развороту МКС. Задача оптимального управления ставится для объекта управления — тела с упругими свойствами. Компоненты управления — импульсные величины, соответствующие моментам сил, создаваемым реальным включением групп РД объединенной двигательной установки МКС по шести каналам управления. В задаче учитываются до 200 упругих мод колебаний конструкции МКС в диапазоне -0...10 Гц. В задачу оптимального управления включено также до 77 фазовых ограничений, соответствующих максимально допустимым значениям упругих нагрузок в 64-х критических интерфейсах АС и РС конструкции МКС. Заметим, что данные по параметрам упругой конструкции МКС, включая коэффициенты нагрузок, предоставлены NASA. Описанная задача оптимального управления была численно решена специалистами ООО «ДАТАДВАНС» (г. Москва) при участии специалистов РКК «Энергия». Решение задачи получено с помощью алгоритмов и программного обеспечения разработки фирмы ООО «ДАТ-АДВАНС». Результатом решения являются как оптимальная последовательность импульсов РД (оптимальная циклограмма РД), так и оптимальная траектория углового движения МКС (Q6(t), <5(t)}, где co(t) вектор угловой скорости МКС относительно ИСК, спроектированный на оси XYZCM базовой СК СМ. Решения получены для шести различных типов разворотов и при этом — для нескольких значений длительностей разворотов Т .
1 1 max
Далее в тексте российская постановка задачи ОРМ и ее численные решения именуются «ОРМ РС». Заметим, что на данном этапе задача синтеза ОРМ РС решена без учета влияния набегающего потока атмосферы на угловое движение МКС.
Выберем следующую последовательность разворотов от осей ОСК к текущим осям XYZCM следующим образом:
• вокруг оси ОYОСК на угол курса в;
• вокруг оси ОZОСК на угол тангажа а;
• вокруг оси ОХСМ на угол крена у.
Траектория оптимального разворота
ОРМ РС для разворота МКС mXVV-pXVV является совокупностью переворота по крену на 360°, разворота по курсу на 174° и отклонения по тангажу от 0° на ~57° с возвратом в 0°. На рис. 3 показаны графики изменения параметров углового управляемого движения МКС во время разворота. На рис. 3, а приведены, сверху вниз, графики углов у, р и а. На рис. 3, б -графики угловых скоростей МКС относительно ИСК в проекциях на оси XCM, YCM и ZCM. Для сопоставления с решением NASA на этом же рисунке приведены также параметры оптимальной траектории ОРМ NASA. Графики угловых скоростей МКС для ОРМ РС, представленные сплошными линиями, имеют скачки, т. е. разрывы гладкости, в моменты
включений/выключений РД МКС. Этого не наблюдается на соответствующих графиках ОРМ АС, полученных решением задачи с «фиктивным» управлением РД. Заметим, что графики ОРМ РС являются результатом численного моделирования управляемого движения МКС, выполненного на стенде РКК «Энергия» путем выдачи фиксированной оптимальной последовательности (циклограммы) импульсов РД, полученной решением задачи ОРМ РС фирмой ООО «ДАТАДВАНС».
Представленная на рис. 3 траектория ОРМ РС получена решением задачи ОРМ для сборки МКС, изображенной на рис. 1.
3. Два подхода к расчету 40-базы АДХ
Как было указано в разд. 1 статьи, наиболее трудный случай при выборе области изменения параметров 4_0-пространства {аа, Ра, ус, 0С} возникает при синтезе оптимальных разворотов МКС на большие углы относительно ОСК, когда оптимальную траекторию разворота предсказать очень трудно.
Рис. 3. Углы оптимального разворота МКС относительно ОСК (а) и угловые скорости относительно ИСК (б):
- - - — оптимальная траектория разворота ОРМ NASA; — оптимальная траектория ОРМ РС, рассчитанная российской стороной без учета аэродинамики
Это видно из описания траектории оптимального разворота в разд. 2. В этой ситуации, чтобы провести синтез оптимального управления с учетом аэродинамики, требуется предварительно рассчитать 4Б-базу АДХ Б4Б во всем диапазоне изменения четырех углов, например, во всех 419 904 узлах КЭ-разбиения Е10.
Для расчетов такого рода РКК «Энергия» применяет два подхода: упрощенный метод и метод расчета с учетом реальной геометрии и затенения. Как указано во введении к статье, геометрия МКС и, следовательно, содержимое базы Б4Б зависят от закона управления 16 приводами вращения панелей СБ и радиаторов во время выполнения разворота. Исследуемая конфигурация МКС представлена на рис. 1.
Во время разворота приводы «Альфа» СБ АС фиксировались в положениях 100 и 255° для правого и левого бортов, соответственно. Приводы «Бета» всех шести панелей СБ АС отслеживали Солнце. Привод «Бета» для СБ АС «1В» фиксировался в положении 230°, привод СБ АС «3А» — в положении 33°. Приводы больших радиаторов АС фиксировались в положениях 75 и 75° для правого и левого бортов. СБ СМ отслеживали Солнце.
Для описанного выше закона вращения/ фиксации поворотных элементов МКС была выполнена процедура расчета узловых АДХ для КЭ-разбиения Е10, описанная в разд. 1 статьи:
• рассчитано множество данных К2{уС , 0'с, Б}, состоящих из КЬ = 703 строк, содержащих все возможные наборы Б поворотов 16 приводов для всех узлов Е10;
• для каждого набора Б(уС , 0'с) рассчитаны 6-векторы АДХ для каждого из !•] = 703 вариантов направления вектора aaer набегающего потока, т. е. 4Б-база Б4Б;
• расчет базы проводился двумя методами:
- упрощенным методом, результат — база
- методом с учетом реальной геометрии и затенения, результат — база БрБ.
4. упрощенный метод расчета 4.0-базы Адх
В основу алгоритма положена расчетная схема, реализующая так называемую зеркально-диффузную аэродинамическую модель обтекания свободномолекулярным потоком набегающих частиц. Суть этой модели заключается в том, что часть молекул набегающего потока претерпевает абсолютно упругое соударение, а остальная
часть молекул в набегающем потоке претерпевает абсолютно неупругое соударение с корпусом МКС [8, с. 12-15]. Разбиение потока на эти две части регулируется коэффициентом диффузности с. В рассматриваемом алгоритме коэффициент диф-фузности с один и тот же для всех элементов конструкции, и его значение по умолчанию равно 0,65, т. е. 65% частиц в потоке претерпевают абсолютно неупругое соударение с корпусом станции, а 35% испытывают абсолютно упругий удар. Коэффициент диффузности был выбран равным 65% по соглашению с NASA.
Для упрощения вычисления аэродинамических коэффициентов каждый элемент конструкции аппроксимируется тремя плоскостями, параллельными плоскостям Oxy, Oxz и Oyz, с заданной площадью и координатами центров давления. Описанные исходные данные были получены из проектных документов МКС от американских коллег.
В соответствии с моделью зеркально-диффузного отражения векторы аэродинамических коэффициентов C*p и Cj для элементарной панели i (2) рассчитываются по следующим формулам:
C = - S
1
{2с A.|(n., а )| а +
1 г IV гJ aer I aer
хар
(3)
+ 4(1 - |(П., Aaer)|(n. , Aaer)n.};
{2cA.|(n. , а )|[r. x а ] +
1 г iv . J aer'' L . aer->
C* =__
T S L 1 1 'V *' ~aer' i l . "aer-i
хар хар
+ 4(1 - c)A.|(n., а )|(n., а )[r. x n.]},
г . aer . aer . .
(4)
где с — коэффициент диффузности; aaer -единичный вектор в направлении против набегающего потока; A. — площадь i-ой аппроксимирующей плоскости; r. — эффективный радиус-вектор из центра масс к центру давления i-ой аппроксимирующей плоскости; n. — единичный вектор нормали к i-ой аппроксимирующей плоскости; S — характерная площадь МКС; £хар — характерный размер МКС. Моментные характеристики определены относительно положения ЦМ.
Описанную выше схему вычисления аэродинамических характеристик широко применяют специалисты NASA для моделирования движения Международной космической станции вокруг ЦМ для верификации алгоритмов управления и синтеза ОРМNASA [6, 8].
На языке С++ был создан модуль расчета АДХ на основе модели аэродинамики, описанной выше, который в качестве исходных данных использует ансамбль строк K2{va 9С, D}. Расчет не учитывает взаимного геометрического затенения элементов конструкции МКС и эффектов интерференции. В модели было использовано 23 различных элемента МКС, т. е. 69 аппроксимирующих плоскостей. Это упрощение позволило рассчитать узловые АДХ во всех 494 209 узлах 4D-разбиения Е10 для МКС с изменяющейся геометрией, изображенной на рис. 1. Расчет АДХ во всех 494 209 узлах занимает не более двух минут на процессоре Intel Core i5 [email protected] GHz. После нормировки и упаковки рассчитанных АДХ в двухбайтные слова соответствующая 4D-база заняла объем 9 084 672 байта, что позволило целиком содержать ее в ОЗУ ЭВМ в алгоритмах моделирования движения и синтеза оптимального управления.
5. Метод расчета АДХ с учетом реальной геометрии и затенения и построения 40-базы для КЭ-разбиения Н10
Более точные результаты расчета аэродинамики могут быть получены по методике, реализованной в программном комплексе Rusat [9], разработанном ИТПМ СО РАН и широко применяемом РКК «Энергия» для расчета АДХ на орбитальном участке полета. При этом учитываются особенности реальной геометрии МКС, эффекты затенения элементов поверхности и интерференции, используется современная модель взаимодействия молекул набегающего потока с различными типами обтекаемых поверхностей.
Однако технология расчета, используемая в программном комплексе Rusat, не позволяет проводить расчеты большого объема при постоянном изменении геометрии объекта. Это связано с тем, что на подготовительном этапе все варианты геометрии сохраняются на жестком диске компьютера и требуют значительных затрат времени, на порядки превышающих время собственно расчета.
В связи с этим был создан новый программный комплекс, выполняющий в т. ч. расчеты свободномолекулярной аэродинамики методом интегрирования по поверхности локальных аэродинамических сил и моментов с учетом эффектов затенения. Для ускорения проверки видимости элементов используется просмотр ячеек 3D-сетки вдоль луча, исходящего из центра
треугольной панели и направленного параллельно набегающему потоку в обратную сторону. Для хранения списка треугольников, пересекающих ячейки сетки, применяется компактная схема в виде массива, упорядоченного по номеру ячейки. Такой способ позволяет обрабатывать сетки с большим числом элементов. Это решает проблему с увеличением объема памяти ЭВМ пропорционально площади поверхности обтекаемого тела и с логарифмической сложностью алгоритма построения списка треугольников, принадлежащих данной ячейке 3^-сетки. Эффекты интерференции при таком подходе не учитываются, но для рассматриваемых конфигураций МКС их вклад не превышает 5%, а их учет требует больших затрат машинного времени.
Этот программный комплекс интегрирован в свободно распространяемую платформу Salome [10], которая в современной графической среде позволяет создавать упрощенную геометрическую 3-О-модель МКС, строить поверхностную расчетную сетку, назначать свойства материалов на сеточной модели, задавать необходимые начальные данные и запускать расчет задачи.
Применяемая схема расчета не требует генерации полной геометрической модели. Проверка затенения выполняется в локальной системе координат каждого неподвижного или трансформируемого сеточного блока. Связи между локальными и глобальной системами координат заданы структурой геометрической модели и параметрами ориентации поворотных элементов МКС, входящих в набор исходных данных.
В рамках используемой методики необходимо знать только локальные коэффициенты обмена импульсом в каждой точке обтекаемой поверхности, зависящие от местного угла атаки. В расчетах принята четырехпараметрическая схема, основанная на модели взаимодействия набегающего потока с поверхностью (модель Ночиллы) [11].
Внешние поверхности МКС моделируются следующими основными типами материалов:
• поверхности модулей — стеклоткань, металл;
• радиаторы — белая эмаль;
• солнечные батареи — стекло.
Температура поверхности МКС принимается равной 300 К. Для высоты орбиты МКС ~400 км в расчетах используются следующие параметры: скорость набегающего потока 7 500 м/с; температура 1 000 K; средний молекулярный вес 16. Угол атаки, скольжения и ориентацию поворотных
элементов конструкции получают из предварительно рассчитанного массива данных К2^С, 6>С, D} для КЭ-разбиения Е10.
Расчет аэродинамики может выполняться как в параллельном режиме для ускорения единичного расчета, так и на кластере путем одновременного запуска большого количества задач на одном ядре через систему очередей.
Упрощенная геометрия сборки МКС, изображенной на рис. 1, состоит из неподвижной части и 13 вращающихся элементов. Общее число панелей модели — 937 134. Среднее время расчета одного варианта — 2,7 с на одном ядре кластерного процессора AMD Opteron(tm) 6378. Общий счет задания из 494 209 вариантов на 64 ядрах потребовал 51 ч машинного времени кластера.
6. вычисление аэродинамических сил и моментов в ходе моделирования оптимального разворота мкС
Для иллюстрации процесса вычисления аэродинамических сил, моментов сил и влияния атмосферы Земли на процесс разворота МКС приведем результаты верификации оптимального разворота, представленного на рис. 3 и для конфигурации МКС, представленной на рис. 1. Верификация
проводилась с помощью моделирующего комплекса РКК «Энергия» МКС-МА («МКС-Модули Автономные»), разработанного одним из авторов данной статьи.
Напомним, что в качестве первого приближения оптимальная циклограмма работы РД МКС и ее оптимальная траектория разворота были вычислены без учета аэродинамики. С помощью моделирования, учитывающего влияние атмосферы Земли, была вычислена траектория разворота МКС продолжительностью 5 500 с при работе РД МКС по оптимальной циклограмме. На рис. 4 представлены графики изменения параметров, входящих в формулы (1) и (2) вычисления аэродинамических моментов сил. Графики получены для следующих начальных данных моделирования: угол Р5 между направлением на Солнце и плоскостью орбиты равен 47,86°; угловое расстояние а5 до орбитального полдня 64,87°; высота орбиты 381,46 км; эксцентриситет орбиты 0,00064; индекс солнечной активности _Р107 = 188,72; индекс геомагнитной возмущенности ар = 24,62. На рис. 5 представлены параметры, требуемые для вычисления АДХ на интервале моделирования, и результат — аэродинамические моменты сил, вычисленные в ходе моделирования 4Б-интерполяцией по «точной» Бе4В и «упрощенной» БББ 4Б-базам.
Рис. 4. Параметры скоростного напора во время разворота МКС: а — плотность атмосферы кг/м3; б — скорость центра масс МКС относительно атмосферы Уд((), км/с; в — скоростной напор Qae, кг/(м-с2)
Рис. 5. Углы-параметры траектории в 40-пространстве во время разворота МКС, а также компоненты
аэродинамического момента сил Та относительно центра масс МКС, н-м, вычисленные 40-интерполяцией: а — углы Vс, 9С (точечная и штрихпунктирная линии, соответственно); б — угол атаки аа и угол скольжения ра (сплошная и пунктирная линии, соответственно); в, г, д — компоненты ТX, Т,, Т,, соответственно, н-м (сплошные линии — моменты сил ТаР, вычисленные интерполяцией по «точной» базе Брв; пунктирные — моменты Та5, вычисленные интерполяцией по «упрощенной» базе Вб4В)
Вертикальные линии на рис. 5, а, б обозначают моменты времени, в которые 4_0-траектория Ъ,(Ь) переходит через границу очередного из 101 4_0-куба разбиения Е10. Также обращает на себя внимание поведение угла атаки аа. В конце разворота должно быть аа = 0°. Значительный уход от нуля, начиная с момента Ь ~ 4 260 с, обусловлен деформацией оптимальной траектории разворота под влиянием аэродинамики, о чем подробнее сказано в разд. 8 статьи.
На рис. 6 представлены АДХ сил и моментов сил относительно базовой СК СМ с началом в центре плоскости стыка узла АО СМ, вычисленные 4^0-интерполяцией по двум 4Л-базам, БР4В и Б^в.
Рис. 5 и 6 иллюстрируют как схожесть в целом АДХ, рассчитанных по точной и «упрощенной» 4.0-базам, так и значительные расхождения, в частности, для АДХ С£, С , СХ и СТ в середине разворота и в окрестности окончания разворота. Это происходит из-за учета в точном расчете взаимных затенений элементов МКС.
Рис. 6. АДХ сил и моментов, вычисленные 40-интерполяцией во время разворота МКС: а, б, в — АДХ сил С X, С I,
соответственно; г, д, е — АДХ моментов сил относительно начала системы координат служебного модуля СX, СI, С?, соответственно
С ?
Примечание. Сплошные линии — АДХ, вычисленные интерполяцией по «точной» базе Бпунктирные линии вычисленные интерполяцией по «упрощенной» базе Б Здесь все АДХ безразмерны.
АДХ,
7. верификация и оценка качества процесса 40-интерполяции Адх
Для оценки качества реализации всего цикла вычисления текущих АДХ МКС с помощью 4Б-интерполяции был проведен
численный эксперимент с использованием данных моделирования оптимального разворота, изображенного на рис. 3. В ходе моделирования вычислялись и записывались с шагом 2 с текущие значения углов поворотов всех приводов СБ и радиаторов
МКС, т. е. ансамбль К2{ус, 9С, Б) строк БЦ.) 16 положений й. всех приводов «вдоль траектории разворота МКС». На рис. 7 приведены графики изменения во времени восьми приводов «Бета» СБ АС.
Используя приведенные данные в качестве исходных, модулем упрощенного расчета АДХ были рассчитаны АДХ МКС через каждые 2 с, т. е. АДХ «вдоль траектории разворота». Анализ показывает, что расхождение между АДХ, рассчитанными «вдоль траектории» по текущим положениям поворотных элементов МКС, и АДХ, рассчитанными 4Б-интерполяцией по текущим положениям векторов набегающего потока атмосферы и направления на Солнце, весьма незначительны.
Для количественной оценки различия двух способов введем следующие обозначения. Пусть С' — максимальное значение
^ ^шах
модуля «силовых» АДХ по оси I, рассчитанное «прямым» способом: С'Р = шах|С'(0|,
0 < Ь < Гшах, I е {X, I, ?}. Аналогично, пусть С^шах — максимальное значение модуля АДХ моментов вокруг оси I, рассчитанное «прямым» способом: СтТшах = шах|СТ.(Ь)|, 0 < Ь < Тшах, I е {X, I, 2). Для оценок расхождения \С'Р(Ь) между текущими значениями С() АДХ, вычисленными «прямым» способом, и значениями СР4В(Ь) АДХ, вычисленными 4Б-интерполяцией, выберем соотношения:
Д£(0 =
Ш) - а
э(Щ
а
•100%
(5)
для силовых АДХ и
дс*(О =
\С1т(Т) С'т 4Д(Т)|
С
•100%
(6)
Ттах
для АДХ моментов.
Как следует из определения, приоритет отдается АДХ, рассчитанным «прямым» способом «вдоль траектории».
Рис. 7. Текущие положения приводов «Бета» СБ АС во время разворота: а — углы СБ АС правого борта МКС (1А — сплошная и ЗА — пунктирная линии); б — углы СБ АС правого борта МКС (1В — сплошная и ЗВ — пунктирная линии); в — углы СБ АС левого борта МКС (2В — сплошная и 4В — пунктирная линии); г — углы СБ АС левого борта МКС (2А — сплошная и 4А — пунктирная линии). Значения углов даны в градусах
Рис. 8. Рассчитанная по соотношениям (3), (4) разница между АДХ, вычисленными вдоль траектории, и АДХ, вычисленными 40-интерполяцией: а, б, в — текущие значения параметров АСХ, АС^, АС^, соответственно, %; г, д, е — текущие значения параметров АСХ, АСТ, АС^, соответственно, %
В нашем примере максимумы рассчитанных «вдоль траектории» АДХ таковы: CX = 427,49; CI = 333,26; CI = 518,10;
Fmax ' 7 Fmax ' 7 Fmax ' 7
CX = 122,60; CT = 1 222,20; CT = 759,19.
Tmax 7 7 Tmax 7 7 Tmax 7
Результаты сравнительного анализа,
проведенного с использованием соотношений (3), (4), приведены на рис. 8.
Максимальные ошибки сравнения составили: АСХ = 3%; ДСГ = 3,4%; & = 1,0%;
гшах 7 гшах ' 7 гшах
ДСХшах = 13,5%; ДСГшах = 1,1%; ДСТшах = 4,5%.
Рис. 9 иллюстрирует локальный характер различий между АДХ, вычисленными «прямым» способом, и 4Б-интерполяцией. Приводятся в увеличенном масштабе графики АДХ (сверху вниз): СХ С^ С2Т и СХ в местах «локальных пиков» графиков рис. 8.
Аналогичным образом, используя исходные данные рис. 7, были рассчитаны АДХ «вдоль траектории разворота» также и «точным» способом. Так же, как и в случае упрощенного способа вычисления АДХ, проведем количественную оценку различия «прямого» расчета АДХ и 4Б-интерполяции по формулам (5) и (6). В этом случае максимумы рассчитанных «вдоль траектории» АДХ таковы: СХ = 397,75; С! = 251,17;
^ гшах ' 7 гшах
& = 337,02; CX = 234,05; CI = 882,21;
Fmax 7 7 Tmax 7 7 Tmax 7 7
CZ = 636,83. Вычисленные максимальные
Tmax
ошибки сравнения с 4D-интерполяцией по базе BP4D составили: ACX = 2,7%;
ACLx = 3,7%; ACLx = 3,2%; A^Lx = 9,1%;
ACTmax = 2,5%; ACTmax = 3,5%.
Проведенный сравнительный анализ и графики рис. 6 и 8 показывают: процедура 4D-интерполяции верно вычисляет АДХ, используя рассчитанную заранее 4D-базу АДХ. При этом ошибки, по сравнению с «прямым» способом вычисления, не превышают 13,5% для «упрощенного» метода и 9,1% — для точного метода расчета АДХ, т. е. лежат в пределах систематических ошибок вычисления АДХ обоими методами.
Рис. 9. Сравнение графиков АДХ в окрестности моментов времени, соответствующих максимальным различиям:
сплошные линии — АДХ, вычисленные 4Б-интерполяцией по базе ББ4В; пунктирные — АДХ, вычисленные «вдоль траектории» упрощенным алгоритмом. Все параметры безразмерны
8. Влияние аэродинамики на траекторию разворота МКС шХУУ-рХУУ
Как видно из рис. 5, величины компонент аэродинамического момента невелики во время разворота. Тем не менее, они оказывают сильное влияние на оптимальную траекторию (см. рис. 3). Для оценки ошибок, которые вносит влияние аэродинамики на номинальную траекторию, выбраны два
скалярных параметра: «ошибка по углу» аегг и «ошибка по скорости» шегг. Пусть требуемые конечные параметры разворота МКС характеризуются кватернионом Орь и угловой скоростью юр а текущие ориентация и угловая скорость — 0ь(0 и ю(Ь), соответственно. Текущая «угловая ошибка» аегг(Ь) — это угол конечного поворота кватерниона Оегг = (ОЬ)* ° Оь(Ь), а текущая «относительная ошибка по скорости»
®от(0 = | - Ю(Ь) | /юорб. ^десь (0РЬ) тернион сопряженный 0Ь, а ш
орб
ква-модуль
угловой скорости движения МКС по околоземной орбите. На рис. 10 представлены графики поведения ошибок а (Ь) и ш (Ь) для номинальной траектории разворота ОРМ РС без учета аэродинамики, а также для траекторий ОРМ РС, возмущенных аэродинамикой, вычисленной двумя способами.
Рис. 10. Поведение ошибок по углу (а) и скорости (б) для трех вариантов моделирования разворота МКС шХУУ-рХУУ:
сплошные линии — моделирование без учета аэродинамики; пунктирные — учет аэродинамики интерполяцией по «упрощенной» базе ; штрихпунктирные — учет аэродинамики интерполяцией по «точной» базе БР4В
Учет аэродинамики сильно деформировал оптимальную траекторию. В конце разворота, в момент ТР = 5 500 с, угол рассогласования возмущенной траектории от цели составил более 120°. Следует заметить, что в момент Ь = 4 402 с угол рассогласования возмущенной траектории от цели был менее 15°, а разница между текущей и целевой угловыми скоростями менее 0,04 °/с. Естественно предположить, что в этот момент времени следует перейти к отслеживанию оптимальной траектории разворота МКС.
Приведенный пример демонстрирует необходимость учета влияния аэродинамики МКС на стадии синтеза оптимального управления, а также необходимость введения дополнительного управления с обратной связью для коррекции неучтенных внешних возмущений.
Заключение
В данной работе представлены начальные результаты исследований, проводимых
РКК «Энергия», по учету сил и моментов, действующих на большую космическую конструкцию, находящуюся на низкой околоземной орбите с высотой ~400 км. Предполагается, что космическая конструкция обладает многочисленными поворотными элементами с разнонаправленными осями вращения. Примерами таких космических объектов являются существующая МКС, а также космические станции будущего.
Предложен метод расчета текущих АДХ и моментов сил аэродинамики при произвольном положении больших космических конструкций относительно набегающего потока и направления на Солнце. Метод разрабатывается для быстрого расчета АДХ при численном синтезе оптимального управления разворотами МКС и при верификации длительных динамических режимов МКС. Представлены различные этапы процесса реализации этого метода на примере моделирования оптимального разворота МКС на большой угол.
Также в процессе моделирования разворота проведено сравнение двух способов вычисления АДХ: упрощенного метода расчета, широко употребляемого коллегами из США, и более «точного» метода с учетом реальной геометрии и затенения, применяемого специалистами РКК «Энергия». При планировании длительных динамических режимов с редкими включениями РД предпочесть следует «точный» метод, несмотря на значительные затраты вычислительных ресурсов.
Показано, что на этапе численного синтеза оптимального управления движением МКС необходимо учитывать влияние атмосферы Земли.
Список литературы
1. Микрин Е.А. Перспективы развития отечественной пилотируемой космонавтики (к 110-летию со дня рождения С.П. Королёва) // Космическая техника и технологии. 2017. № 1(16). С. 5-11.
2. ГОСТ 20058-80. Динамика летательных аппаратов в атмосфере. Термины, определения и обозначения. М.: «Издательство стандартов», 1981.
3. Микеладзе В.Г., Титов В.И. Основные геометрические и аэродинамические характеристики самолетов и ракет. Справочник. М.: Изд-во «Машиностроение», 1990. 149 с.
4. ГОСТ Р 25645.166-2004. Атмосфера Земли верхняя. М.: ИПК «Издательство стандартов», 2004.
5. Зенкевич О. Метод конечных элементов в технике. Пер. с английского. М.: «Мир», 1975. 543 с.
6. Bedrossian N, Bhatt S, Kang W., Ross M. Zero-propellant maneuver guidance // IEEE Control Systems magazine. October 2009. V. 29. № 5. P. 53-73.
7. Bhatt S., Bedrossian N., Nguyen L. Optimal propellant maneuver flight demonstrations on ISS // AIAA Guidance, Navigation, and Control Conference 2013. 19-22 August 2013, Boston, Massachusetts, USA (AIAA 2013-5027). Режим доступа: http://dx.doi.org/10.2514/62013-5027 (дата обращения 06.09.2017 г.).
8. Bhatt S. Optimal reorientation of spacecraft using only control moment gyroscopes // Master's Thesis, Rice University, Houston, TX, May 2007. P. 110. Режим доступа:
http://www.caam.rice.edu/tech_reports/2007/ TR07-08.pdf (дата обращения 06.09.2017 г.).
9. Кашковский А.В, Ващенков П.В., Иванов М.С. Программная система для расчета аэродинамики космических аппаратов // Теплофизика и аэромеханика. 2008. Т. 15. № 1. С. 79-91.
10. The Open Source Integration Platform for Numerical Simulation. Copyright © 2005-2017 OPEN CASCADE. Режим доступа: http:// www.salome-platform.org/ (дата обращения 06.09.2017 г.).
11. Freedlender O.G., Nikiforov A.P. Modelling aerodynamic atmospheric effects on the space vehicle based on test data // Proceedings of 2nd Int. Symp. on Environment testing for space programmes, 1993. P. 307-312. Статья поступила в редакцию 23.08.2017 г.
Reference
1. Mikrin E.A. Perspektivy razvitiya otechestvennoi pilotiruemoi kosmonavtiki (k 110-letiyu so dnya rozhdeniya S.P. Koroleva) [Outlook for our country's manned spaceflight development (to mark the 110th anniversary of S.P. Korolev)]. Kosmicheskaya tekhnika i tekhnologii, 2017, no. 1(16), pp. 5-11.
2. GOST 20058-80. Dinamika letatel'nykh apparatov v atmosfere. Terminy, opredeleniya i oboznacheniya [Dynamics of flying vehicles in the atmosphere. Terms, definitions and acronyms]. Moscow, Izdatel'stvo standartovpulb., 1981.
3. Mikeladze V.G., Titov V.I. Osnovnye geometricheskie i aerodinamicheskie kharakteristiki samoletov i raket. Spravochnik [Basic geometric and aerodynamic properties of airplanes and rockets]. Moscow, Mashinostroeniepubl., 1990. 149p.
4. GOST R 25645.166-2004. Atmosfera Zemli verkhnyaya [Upper atmosphere of Earth]. Moscow, IPK Izdatel'stvo standartov publ., 2004.
5. Zenkevich O. Metod konechnykh elementov v tekhnike [Finite elements method in engineering]. Moscow, Mir publ., 1975. 543 p.
6. Bedrossian N, Bhatt S, Kang W, Ross M. Zero-propellant maneuver guidance. IEEE Control Systems magazine, October 2009, vol. 29, no. 5, pp. 53-73.
7. Bhatt S, Bedrossian N, Nguyen L. Optimal propellant maneuver flight demonstrations on ISS. AIAA Guidance, Navigation, and Control Conference 2013, 19-22 August 2013, Boston, Massachusetts, USA (AIAA 2013-5027). Available at: http://dx.doi.org/10.2514/62013-5027 (accessed 06.09.2017).
8. Bhatt S. Optimal reorientation of spacecraft using only control moment gyroscopes. Master's Thesis, Rice University, Houston, TX, May 2007, p. 110. Available at: http: //www.caam.rice.edu/tech_reports/2007/TR07-08.pdf (accessed 06.09.2017).
9. Kashkovskii A.V, Vashchenkov P.V., Ivanov M.S. Programmnaya sistema dlya rascheta aerodinamiki kosmicheskikh apparatov [Software system for computing aerodynamics of spacecraft]. Teplofizika i aeromekhanika, 2008, vol. 15, no. 1, pp. 79-91.
10. The open source integration platform for numerical simulation. Copyright © 2005-2017 OPEN CASCADE. Available at: http://www.salome-platform.org/ (accessed 06.09.2017).
11. Freedlender O.G., Nikiforov A.P. Modelling aerodynamic atmospheric effects on the space vehicle based on test data. Proceedings of 2nd International Symposium on Environment testing for space programmes, 1993. P. 307-312.
