Применение суперкомпьютерных технологий для исследования авиационных катастроф

С учетом исключительных знаний и компетен­ций, которые есть у Концерна ВКО «Алмаз - Антей» и его предприятий, и понимания необ­ходимости точного и всестороннего исследо­вания авиационной катастрофы малазийского самолета «Боинг» (рейс МН-17), произошед­шей в июле 2014 г. над территорией Украины, руководством Концерна в мае 2015 г. были поставлены следующие задачи:

• создать точную математическую 3.0-мо­дель для анализа деталей летного происшествия;
• по результатам моделирования полу­чить данные, необходимые для проведения натурных экспериментов;
• определить возможные точки старта зенитной управляемой ракеты (ЗУР).

Для выполнения поставленных задач, требующих исключительной вычислитель­ной точности, был разработан программный комплекс моделирования авиационного про­исшествия с применением суперкомпьютера «Орфей-К» (рис. 1) [1], состоящий из пяти программных модулей (ПМ).

Программный модуль расчета разлета осколков боевой части (БЧ) ЗУР производит расчет на основании годографа, определяю­щего зависимость скорости движения по всем поражающим элементам (легкая и тяжелая фракции, а также элементы корпуса БЧ) от их направления [2].

Модель разлета осколков позволяет вычислить скорость движения и координа­ты каждого из них в любой момент времени. Осколки в программном комплексе моделиру­ют с помощью материальных точек, их движе­ние описывается кинематическим уравнением прямолинейного равномерного движения:

где (x_i, y_i, z_i) - координаты осколков в момент времени t;

(x_0i, y_0i, z_0i) - координаты осколков в на­чальный момент времени t = 0;

(V_xi, V_yi, V_zi) - компоненты скоростей дви­жения осколков.

Скорость разлета каждого из осколков за­дает годограф, описывающий конкретную БЧ.

ПМ спроектирован таким образом, что позволяет работать практически с любой ЗУР, так как считывает данные о годографе из кон­фигурационного файла, сохраняет их в специализированной базе данных и использует для расчета в модели.

Программный модуль визуализации ре­зультатов моделирования предназначен для трехмерного отображения результатов моде­лирования: поверхностей корпусов самолетов («Боинг-777-200 ER» и Ил-86), поражения ча­стей конструкции самолета (кабины пилота, остекления кабины, крыльев, двигателей, хво­стового оперения, пассажирского салона и т. д.) различными фракциями поражающих элемен­тов, углов входа поражающих элементов в ча­сти конструкции самолета, направления дви­жения поражающих элементов, точек старта ЗУР с привязкой к местности и др. (рис. 2).

Программный модуль «динамика-ста­тика» предназначен для выполнения тру­доемких вычислений на суперкомпьютере «Орфей-К» и служит математической основой для натурных исследований. Очевидно, что проведение эксперимента в реальных услови­ях, когда цель и ЗУР подвижны, невозможно. В этом случае ПМ ведет поиск статического положения ЗУР, которое показывает повреж­дения на поверхности, максимально схожие с теми, которые нанесет ЗУР в движении [3].

Смысл алгоритма решения задачи «По­иск статики по динамике» в следующем. В некоторой области пространства производят подрыв БЧ ЗУР в «динамике», когда БЧ и цель движутся. Вычисляют траекторию движения каждого из осколков и точки пересечения тра­екторий с поверхностью самолета. Найденные точки сохраняют в массиве данных, каждый элемент которого представляет собой число пробоин, приходящихся на площадку с но­мером, соответствующим индексу элемента.

Массив, хранящий число пробоин в каждой из площадок самолета для подрыва БЧ в динами­ке, назовем эталоном.

Далее скорость цели и ЗУР обнуляют и «выполняют» подрывы БЧ ЗУР в различных точках пространства с разными углами ази­мута и места. При каждом подрыве опреде­ляют точки пересечения траекторий осколков с поверхностью самолета. Затем формируют массив, хранящий число пробоин в каждой из площадок. Будем называть такие массивы ап­проксимациями (приближениями).

Для всех найденных аппроксимаций вычисляют ошибку Е, количественно харак­теризующую различия между массивом ап­проксимации и эталоном и использующуюся в качестве критерия отбора массива аппроксима­ции, наиболее совпадающего с эталоном. Зна­чение ошибки Е для него - самое низкое среди значений остальных аппроксимаций.

В итоге решение задачи «Поиск статики по динамике» сводят к численному решению задачи оптимизации

где  - пятимерный  вектор;

(x, y, z) - декартовы координаты центра БЧ ЗУР;

а и е — азимут и угол места БЧ ЗУР;

 - искомое решение, доставляющее минимум функционалу .

Достоверность результатов оценивают в соответствии со следующим выражением:

Программный модуль расчета траекто­рий полета ЗУР предназначен для расчета тра­екторий полета ЗУР от точки стояния до цели. Возможные точки старта ПМ можно опреде­лить с помощью суперкомпьютера «Орфей-К».

Программный модуль управления спе­циализированной базой данных предназначен для хранения информации о любых цифровых картах местности, моделях, описывающих цели (в нашем случае самолеты «Боинг-777-200 ER» и Ил-86), моделях любых ЗУР. ПМ также со­бирает и обрабатывает результаты расчетов на суперкомпьютере «Орфей-К», например обла­сти старта ЗУР с привязкой к цифровой карте местности, задачи «динамика-статика». База данных имеет открытую масштабируемую архитектуру и позволяет работать с широким набором исходных данных как по целям (на­пример, F-16, F-22, Rafale, Eurofighter Typhoon, AGM 158 JASSM и др.), так и по ЗУР (48Н6Е, 9М317, 9М96, 9М100 и др.).

Проведение вычислительных экспериментов

Разработка трехмерной модели авиационного происшествия была начата с создания модели разлета осколков БЧ.

Скорость каждого из осколков задали годографом, начальные положения осколков соответствовали координатам на поверхности БЧ, характеристики которой были предоставлены предприятием-изготовителем. На рис. 3 представлены компьютерные модели разлета осколков БЧ ЗУР 9М38 и 9М38М1. На первом этапе была разработана простая модель их раз­лета из «точки» (рис. 3, а, б). Модель [5] пред­полагает, что в начальный момент времени все осколки расположены в одной точке:

Далее разработанная модель была усо­вершенствована - разлет осколков произво­дился с «отрезка» (рис. 3, в, г). Для этой моде­ли осколки в начальный момент времени t = 0 расположены на отрезке АВ с концами (x^A, y^A, z^A) и (x^B, y^B, z^B):

Осколки распределены по отрезку АВ следующим образом. Отрезок разделили на равные части. Число осколков на каждой из них выбирали в соответствии с заданной плот­ностью распределения:

где m_i - плотность распределения осколков;

М - число подынтервалов, составляющих отрезок АВ;

Р - число поворотов векторов скоростей осколков вокруг отрезка AB при повороте на полный угол, 360°.

После предоставления данных о геометрии БЧ была создана окончательная модель разлета осколков с боковой поверхности эллиптического цилиндра, форма которого практически совпада­ет с формой БЧ ЗУР (рис. 3, д, е). В данной мо­дели осколки располагают вдоль отрезка в соот­ветствии с выражениями (4) и (5). Для получения пространственного распределения осколков от­резок вращается вокруг центра цилиндра. В ито­ге были получены все типы разлета осколков - от точечного до цилиндрического.

В ходе численных экспериментов уста­новлено, что необходимо учитывать форму БЧ ЗУР, так как она влияет на характер пора­жения цели из-за близкого расположения БЧ к поверхности цели. Варианты модели разлета осколков (см. рис. 3) отличаются по форме об­лака разлетающихся осколков.

Модель поверхности самолета «Бо­инг-777-200 ER» изначально была получена из открытых источников в формате 3D Max. Но эта модель не соответствовала реальной моде­ли самолета, более того, недопустимое отличие в форме кабины не позволяло провести полно­ценный вычислительный эксперимент. Для про­ведения исследования необходимо было сде­лать интерактивную модель самолета для того, чтобы видеть, какие поражения наносят оскол­ки БЧ различным элементам конструкции самолета. После уточнения экспертами Кон­церна деталей чертежей модель была скоррек­тирована (рис. 4).

Модель, разработанная в 3D Max, сфор­мирована из большого числа примитивных элементов, поэтому для обеспечения прием­лемого быстродействия необходимо иметь модель, состоящую приблизительно из не­скольких десятков тысяч элементарных тре­угольников. Для ее создания полученную ин­формацию выгружают в специальный файл (хранится в специализированной базе данных), описывающий вершины треугольников и пред­назначенный для обработки в программных модулях визуализации средствами OpenGL и «динамика-статика».

Подготовка к первому натурному экс­перименту и определила постановку первой трудоемкой задачи «Поиск статики по дина­мике». Иными словами, нужно найти такое статическое положение БЧ ЗУР («статика»), повреждения на поверхности цели от которого максимально совпадают с повреждениями на поверхности цели в динамике («динамика»), когда ЗУР и цель находятся в движении.

При выполнении численных расчетов для первого этапа натурного эксперимента был выбран следующий критерий для мини­мизации:

где М - число площадок, составляющих по­верхность самолета «Боинг-777-200 ER»;

е - массив, хранящий число пробоин, приходящихся на каждую площадку в «дина­мике»;

а - массив, хранящий число пробоин, при­ходящихся на каждую площадку в «статике».

Число Е в этом критерии соответствует количеству несовпавших повреждений. Выбор такого критерия обусловлен его очевидностью и простотой реализации.

Рассмотрим выражение (6) подробнее. Для самолета «Боинг-777-200 ER» M = 5365 (число площадок, реализованных в ПМ визу­ализации), общее число поражающих элемен­тов 10 763. Таким образом, эталонный вектор (для первой фазы натурного эксперимента за­полненный значениями при подрыве ЗУР от­носительно цели с углами 62° (азимутальный угол) и 22° (угол места)) будет иметь вид:

Общее число векторов, с которыми по формуле (6) сравнивают эталонный век­тор, определено в формуле (7) и составляет 64 800 000 возможных положений ЗУР в про­странстве. Достоверность результатов перебо­ра оценивают выражением (2).

Для расчетов, обеспечивающих прове­дение первого этапа натурного эксперимента (рис. 5), были приняты следующие условия: го­дограф - ЗУР 9М38М1 с БЧ 9Н314М; поверхность самолета - «Боинг-777-200 ER»; углы ЗУР относительно цели - азимутальный угол 62° и угол места 22° (данные параметры соот­ветствуют версии Концерна о событиях авиа­ционной катастрофы).

Для решения подобных задач нужно перебрать множество вариантов возможных положений ЗУР и их сравнений с эталонным значением. При этом варианты рассматривают по повреждениям в «статике», а эталоном счи­тают повреждения в «динамике» [4] (рис. 6).

Положение БЧ в пространстве определя­ют с помощью пяти параметров: трех декарто­вых координат, угла места и азимута (рис. 7).

Метод полного перебора вычисли­тельно трудоемкий и требует значительных временнь1х затрат. Так, для расчета одного взрыва БЧ требуется около 2 с (данный пока­затель получен экспериментальными замера­ми на «Орфей-К» [6]). Перебор на «грубой» сетке (АРМ с процессором Intel Xeon E5620 2,4 ГГц) для расстояния 10 м, с шагом 0,5 м и для углов с шагом 1 м и поворотом на 90° зай­мет около 4 лет и составит приблизительно 60 млн вариантов:

V_{кол-во вариантов перебора} = V_X V_Y V_Z V_A V_M            (7)

где V_X - количество положений в простран­стве по оси Х, в нашем случае 10 м/0,5 = = 20 вариантов;

V_Y - количество положений в простран­стве по оси Y, в нашем случае 10 м/0,5 = = 20 вариантов;

V_Z - количество положений в простран­стве по оси Z, в нашем случае 10 м/0,5 = = 20 вариантов;

V_A - количество положений по азимуту, в нашем случае 90°/1 ° = 90 вариантов;

V_M - количество положений по углу ме­ста, в нашем случае 90°/1 ° = 90 вариантов. Итого:

Почти такое же время необходимо для выполнения перебора на «тонкой» сетке (рас­стояние 1,2 м с шагом 0,1 м и поворот на угол 20° с шагом 0,1 м). Перебрать надо приблизи­тельно 70 млн вариантов (значение получено из формул (7) и (8)). Перебор на «грубой» сетке необходим для определения области простран­ства, в которой зафиксировано аппроксимиру­ющее положение. Затем на «тонкой» сетке в заданной области пространства с требуемой точностью можно вычислить наиболее под­ходящее положение. В итоге расчет занял бы около 8 лет.

Сократить время выполнения полного перебора удалось путем проведения расчетов на суперкомпьютере «Орфей-К» [5]. Для по­строения параллельного алгоритма была использована технология MPI (Message Passing Interface - наиболее распространенный ин­терфейс передачи сообщений при параллель­ном программировании [7]), так как именно такой подход позволяет оперативно создавать параллельные программные комплексы. Вре­мя полного перебора как на «грубой» сетке, так и на «тонкой» сократилось до 12 суток - алгоритм оказался идеальным для распарал­леливания. Возможность распараллеливания расчетов обусловлена тем, что расчеты по каж­дой анализируемой точке пространства можно проводить независимо друг от друга, поэтому дополнительные меры для доведения алго­ритма к параллельной архитектуре не нужны, так как он получился естественным образом

при раскладывании расчетов на доступные в «Орфей-К» вычислительные потоки:

где T_{полною набора} = 129 600 000 (на «Грубой» сетке + 138 240 000 (на «тонкой» сетке) = 267 840 000 с; количество вычислительных ядер (потоков) - 256;

Т_{перенабора супер ЭВМ} = 26 780 000/256 ≈ 12,1 су­ток.

Расчет был выполнен, найденные пара­метры положения БЧ в «статике» позволили провести первый этап натурного эксперимен­та. В качестве эталона был принят вектор, состоящий из номеров площадок, на которые разделили интерактивную модель поверхности самолета, и количества осколков, попавших в данную площадку. Вектор имел значения, полученные при подрыве БЧ ЗУР по версии Концерна (углы ЗУР относительно цели - ази­мутальный угол 62° и угол места 22°). Ком­пьютерная модель показывала, что параметры рассчитанного статического положения ЗУР должны совпасть по результатам на 96,77 %.

Проведенный вычислительный эксперимент подтвердил полученные расчеты (рис. 8).

В дальнейшем комплекс параллельных программ был усовершенствован: трудоемкие вычисления были вынесены в конфигурацион­ные файлы для ПЛИС.

Кроме получения картины распределе­ния осколков на поверхности самолета экс­перимент установил характер проникнове­ния осколков внутрь самолета, особенности деформации осколков, особенности пораже­ния поверхности частицами разной фракции и многое др.

Одним из важнейших результатов перво­го этапа натурного эксперимента стало под­тверждение факта поражения элементов левого двигателя (ориентацию элементов конструкции самолета определили по направлению дви­жения самолета) и левого крыла. При этом рас­пределение повреждений полностью совпало с расчетами на математической модели, что еще раз доказало правильность выводов экспертов Концерна относительно направления запуска ракеты (рис. 9, 10). Следует отметить, что все результаты натурного эксперимента были по­лучены в ходе всего лишь одного подрыва БЧ с параметрами, вычисленными на суперком­пьютере.

После первого этапа натурного экспери­мента было принято решение провести второй этап (см. рис. 10) и использовать реальный самолет в качестве цели. С учетом текущих возможностей был выбран Ил-86. Таким об­разом, условия для расчетов, обеспечива­ющих проведение второго этапа натурного эксперимента, были следующие: годограф - ЗУР 9М38М1 с БЧ 9Н314М; поверхность самолета - Ил-86; углы ЗУР относительно цели - азимутальный угол 17° и угол места 7° (данные параметры соответствуют версии Международной комиссии).

Очевидно, что и второй этап эксперимен­та провести в реальных условиях, когда цель и ЗУР находятся в движении, невозможно. Снова возникла необходимость решения за­дачи «Поиск статики по динамике». Для про­ведения компьютерного эксперимента была создана цифровая модель поверхности само­лета Ил-86 (рис. 11).

При выполнении расчетов для второго этапа натурного эксперимента был выбран критерий

В данном случае была реализована идея (подтвержденная результатами расчетов) о получении распределения числа осколков на поверхности кабины в статике, визуально похожего на распределение осколков в динамике. Действительно, если проанализировать плот­ности распределения осколков, приходящихся на площадку, для статического и динамическо­го случаев, то заметно, что картины распре­деления плотности осколков имеют внешнее сходство (рис. 12).

Компьютерный эксперимент проводи­ли на готовом программном комплексе. Из-за ограниченности во времени на пошагово сгу­щающихся сетках было перебрано прибли­зительно 49 млн вариантов. Наилучший по выбранному критерию вариант был найден в последнем «запуске», при этом число пробо­ин E = 85, что дает достоверность η = 99,2 % (рис. 13).

Результаты компьютерного эксперимен­та были подтверждены. На рис. 13 видно, что картины площадок с числом повреждений бо­лее 20 в точности совпадают для статическо­го (аппроксимация) и динамического (эталон) случаев.

При анализе результата, полученного в ходе моделирования на суперкомпьютере, (рис. 14) видно, что как в «статике», так и в «динамике» большое количество осколков по­падает в исследуемое остекление кабины с ле­вой стороны. Кроме того, углы входа основной массы осколков как в «статике», так и в «дина­мике» практически совпадают.

Как и в первом, так и во втором натур­ном эксперименте была использована ЗУР 9М38М1, особенность которой заключается в применении двутавровых осколков в качестве тяжелой фракции. Натурные эксперименты показали, что эти осколки всегда оставляют на обшивке самолета следы в форме бабочки. Таких пробоин на обломках поверхности само­лета «Боинг-777-200 ER» обнаружено не было.

Особое внимание, как и на первом эта­пе натурного эксперимента, было уделено по­вреждениям левого двигателя и левого кры­ла (рис. 15). Компьютерное моделирование, с помощью которого решалась задача «Поиск статики по динамике», однозначно говорило о том, что повреждений на левом двигателе и крыле при таких углах подхода нет. Важ­но подчеркнуть, что данных повреждений не было ни в одном из 49 млн вариантов пере­бранных положений ЗУР и результатами вто­рого эксперимента это было подтверждено: щит, имитирующий левое крыло, не был задет ни одним из более чем 10 тыс. осколков.

Результаты двух натурных эксперимен­тов полностью подтвердили правильность соз­данной компьютерной модели (рис. 16). Это позволило перейти к решению второй трудоемкой задачи - расчету всех возможных точек старта ЗУР с целью поиска наиболее вероят­ной области старта ЗУР по известным углам азимута и места.

Для определения траектории полета раке­ты и соответственно района ее запуска у пред­приятия-изготовителя был получен комплекс программ моделирования изделий 9М38. Дан­ные модели были разработаны еще в середине 1980-х гг. и реализованы на языке программи­рования «Фортран». В этих моделях учтено мак­симальное количество параметров, влияющих на условия встречи ЗУР с целью - подробные частные модели двигательной установки, ав­топилота, головки самонаведения и др. Учи­тывались эффективная площадь рассеяния и линейные размеры цели, ее скорость и другие па­раметры движения. Модель была откалибрована на огромном количестве натурных испытаний.

Для того, чтобы обеспечить проведение имитационного моделирования в сжатые сро­ки, был модифицирован исходный программ­ный код, а также написана отдельная оболочка для обеспечения запуска параллельных расчетов на суперкомпьютере с использованием тех­нологии MPI. Для выполнения оперативных расчетов также был создан отдельный модуль с графическим интерфейсом, на котором мож­но в реальном времени получать значения для конкретных параметров движения зенитной управляемой ракеты и цели.

Возможность указания точек старта свя­зана с соответствиями между параметрами подлета ЗУР к цели и координатами ее стар­та, следовательно, имея параметры пространственного положения БЧ ЗУР относительно самолета и географические координаты пора­жения цели, всегда можно вычислить область на карте, с которой был выполнен старт зенит­ной управляемой ракеты.

После написания MPI-версии програм­мы за 6 суток был просчитан (с интервалом в 10 м, для ЗУР 9М38) возможный район пуска (его площадь, 4900 кв. км, определена дальней границей зоны поражения). В расчет включе­на стрельба как в переднюю полусферу, так и «вдогон». В ходе расчетов была установлена возможная область старта ЗУР 9М38 (рис. 17).

Вывод

Использование разработанного программно­го комплекса позволило объективно и на на­учной основе провести исследование причин авиакатастрофы малазийского самолета «Бо­инг-777-200 ER». Впервые в истории иссле­дования причин авиационного происшествия была проведена столь масштабная и скрупу­лезная работа: аналитические исследования экспертов АО «Концерн ВКО «Алмаз - Ан­тей», проведение натурных экспериментов, математическое моделирование и расчеты на суперкомпьютере. По результатам ана­лиза сделаны следующие выводы: самолет «Боинг-777-200 ER», следовавший рейсом МН-17, был поражен ракетой 9М38; место пуска находилось в районе, расположенном южнее населенного пункта Зарощенское.