Балістична ракета 600 км+ на С++

💡 Усі статті, обговорення, новини про оборонні технології — в одному місці. Приєднуйтесь до DefTech спільноти!

Вітаю спільноту. Мене звати Олексій, я цікавлюсь балістикою і програмуванням. Пишу в Qt, тому що там легко інтегрувати графічний інтерфейс в потужне ядро обчислень C++ для майже будь якої платформи. В Qt зробив для смартфона систему злиття сигналів сенсорів для визначення орієнтації (IMU fusion) на фільтрі Махоні, а потім його замінив на власний фільтр Калмана (ESKF).

Одним із перших моїх складних проєктів в програмуванні було моделювання польоту мінометного пострілу. Просто з цікавості зробив щось на кшталт балістичного калькулятора і віддав знайомим в ЗСУ. Один з них сказав, що для міномета не потрібен, але снайпери деякі розрахунки досі роблять в Excel — їм потрібно зробити додаток для Android, щоб обчислювати це автоматично. Без подробиць скажу, що це школа підготовки снайперів і жоден з існуючих комерційних додатків не міг надати потрібні опції аби втілити їхні ідеї. Так народилась моя насправді серйозна програма з балістики для снайперської зброї — кастомний балістичний калькулятор CalcTheBullet. В ньому ми зробили дійсно унікальні речі, наприклад, врахували вертикальний градієнт вітру і реалізували режим стрільби з закритої позиції без прямої видимості цілі по наведенню з дрона. На YouTube за цією назвою можна знайти відео з демонстрацією роботи програми при здійсненні таких пострілів на дуже далекій дистанції.

Маючи певний бекграунд я задався ціллю змоделювати політ балістичної ракети, цікаво було «чисельно побачити» як вона стартує, яку швидкість набирає, як нею прицілитись і куди можна влучити. Якщо і вам цікаво, то далі поясню з прикладами.

Балістика на C++ в Qt Framework

В якості прототипа моєї ракети я взяв з Вікіпедії та інших відкритих джерел дані по твердопаливній ракеті Грім:

• діаметр — 0,92 метри;
• повна маса — 3780 кг;
• бойова частина — 500 кг;
• суха маса — 1130 кг;
• питомий імпульс — 264 секунди;
• тяга — 171600 Н;
• час вигоряння палива 40 секунд;

На фото видно, що у неї є хвостові керуючі поверхні, отже в щільних шарах атмосфери можна маневрувати.

Для визначення таких характеристик як моменти інерції, центр мас та розмір керуючих поверхонь зробив 3D креслення в AutoCAD.

Перед початком польоту треба змоделювати основні параметри. Аби наблизитись до моделі руху справжньої ракети треба врахувати динамічну зміну маси через вигоряння палива, розробити модель атмосфери змінної щільності і температури, змоделювати гравітацію та неідеально круглу форму Землі, знайти дані по опору тіла зі схожою аеродинамікою, розробити систему координат для суборбітального польоту, а також метод інтегрування. Я звичайно не одразу до цього прийшов, але ось що вийшло.

Атмосфера. Для висот до 45 км використовується стандартна модель NASA. Вище, де тиск незначний, застосовано шарувату експоненціальну модель атмосфери, а межі шарів узгоджені, щоб забезпечити безперервність функції щільності на висотах 100 км і вище.

Форма планети. Адже Земля не є ідеально круглою, то для точності обчислення відстаней і азимуту використовується топографічна система координат WGS84, що точно описує форму планети і є стандартом для GPS-навігації, картографії, геодезії, використовуючи так званий референц-еліпсоїд, що задається екваторіальним і полярним радіусами.

Гравітаційна модель J2. У точній навігації класичної формули Ньютона для кулі недостатньо. Екваторіальне здуття планети створює надлишок маси на екваторі, який працює як додатковий пояс тяжіння. Через нього супутники на орбіті відчувають нерівномірну силу, яка постійно змінює їх траєкторію. Це моделює динамічний коефіцієнт форми J2.

Системи координат. Вже згадана WGS84 — геодезична для роботи з картою. ECEF (Earth-Centered, Earth-Fixed) — система координат, яка обертається разом із Землею. Це стандарт для навігації і балістики, використовується як основна для всіх обчислень. Локальна зв’язана система координат (Body Frame) для аеродинаміки: кут атаки, підйомна сила, тяга. NED (North East Down) — для локальних обчислень відносно поверхні Землі. Трансформація між ними відбувається за допомогою матриць повороту бібліотеки лінійної алгебри Eigen для C++.

Аеродинаміка. У вільному доступі реальних даних аеродинамічного опору ракет знайти неможливо, тому внесено пресет у вигляді балістичного коефіцієнта стандартного снаряду G7 змінної маси з добре відомими табульованими динамічними коефіцієнтами опору залежно від числа Маха.

Метод чисельного інтегрування. Метод Хойна (RK2) з кроком 0,001 с, що забезпечує стабільність розрахунків на високих швидкостях. Компроміс між швидкістю і точністю. Швидкість виконання коду з Eigen у мене чомусь дуже повільна в debug режимі, але добре оптимізується в режимі release.

Архітектура.

Головний диригент, що спілкується з графічним інтерфейсом, передає на обчислення, приймає сигнали і керує класами — AppController.

В підпорядкуванні у нього два головних високонавантажених обчислювача фізичних моделей і стану системи:
Class Calculations. Служить для попереднього обрахунку матриці чутливості запуску. Стартує в окремому потоці і сканує всі траєкторії з можливих кутів запуску — виводить дані на екран у вигляді таблиці: для кожної відстань, час польоту, швидкість біля цілі, апогей тощо. Загалом 91 траєкторія.

Class SimulationWorker. Окремий потік і симуляція польоту в реальному часі. Обробляє цикл while(altitude>0), всередині два таймери: перший QElapsedTimer, який рахує фізичний час польоту, другий таймер інтегратора, який доганяє перший і узгоджує швидкість виконання коду. Якщо є зазор, то потік іде спати, аби не навантажувати даремно процесор. Цикл не обертається швидше ніж реальний таймер, а виконання може сильно запізнюватись для можливості покрокового дебагу. Зроблено це з наступною метою.

Class ImpactEstimator. Коли ракета проходить апогей (найвищу точку траєкторії) вона має знати куди треба здійснити маневр, щоб влучити в задану ціль. Це система передбачення місця падіння. З частотою 20 Герц від поточної точки і до землі йде розрахунок похибки влучання, передається сигнал чи маневр потрібен.

Class CorrectionManager. Вмикається, коли входимо в щільні шари атмосфери. Отримує сигнал про похибку влучання і здійснює керування хвостовими поверхнями. Це призводить до повороту тіла ракети проти потоку повітря, що збільшує опір, народжує підйомну силу, що викривляє траєкторію. Чим нижче спускається ракета, тим вище динамічний тиск і ефективніше викривлення. Таким чином відбувається маневрування балістичної ракети, що заважає ППО передбачити точку перехоплення.

Перелічу інші системи.

GeoMath — математична модель Землі. Тут відбувається трансформація координат одна в іншу, рахуються геодезичні відстані і кути, моделюється гравітація.

MissileState, MissilePhysicsModel — аеродинаміка, інтегрування (розрахунок швидкості, прискорення, поточного положення, орієнтації), сил опору, коріоліса, тяги, відцентрової, оновлення матриць повороту тощо.

TrajectoryRecorder — видає телеметрію польоту в CSV файл і консоль розробника, також хочу ним виводити на екран дані польоту в реальному часі, які в UI ще не видно.

Class Calculations і SimulationWorker використовують одні й ті самі модулі, але в різних конвеєрах. Якщо взяти одну і ту саму траєкторію і прогнати її через ці два різні алгоритми, то різниця в результатах складає ~0,008%. Ніби небагато, але в абсолютних числах на 600 км це аж 50 метрів.

Розглянемо приклад обчислень та як це виглядає.

Мапа вибору цілі

Програма запускається у вигляді інтерактивної мапи, де обирається точка запуска — окреслюється коло з радіусом 670 км, що відповідає максимальній дистанції польоту внесеного прототипу на захід. Через ефект Етвеша дальність польоту на схід більше.

Обирається точка влучання, автоматично розраховується відстань та азимут на ціль. Розрахунок відстані відповідає системі WGS84 і завдяки формулі Вінсенті дає похибку визначення відстані на поверхні Землі 0,5 мм.

Матриця чутливості

Після натискання кнопки «Підготовка» відбувається обрахунок можливих траєкторій польоту, які ми отримаємо при відхиленні вектора тяги двигуна під час старту від 6 до 15 градусів з кроком 0,1 градус. Це займає приблизно 30 секунд на стандартному ПК.

Алгоритм маневру запуска: вертикальний старт протягом двох секунд з наступним відхиленням вектору тяги на заданий кут по азимуту цілі. В цьому вікні представлена матриця чутливості траєкторій польоту в напрямку заданої на попередньому екрані цілі (загалом 91 штука).

Представлені наступні дані.

Тяга — відхилення вектору тяги двигуна, градуси.

Відстань — дистанція балістичного польоту в некерованому режимі, метри.

Апогей — максимальна висота від поверхні еліпсоїда, метри.

Азимут — горизонтальний кут заходу на точку падіння, радіан.

Тангаж — вертикальний кут падіння, радіан.

Швидкість — швидкість біля цілі, м/с

Час — час польоту, секунди.

Корекція сили Коріоліса

Наступне вікно пропонує скоригувати азимут запуску і зменшити бокове відхилення, викликане силою Коріоліса, яке в цьому випадку складає майже 11 км. Якщо ввести запропонований азимут, то відхилення зменшиться до метрів.

Після погодження корекції азимута запускається симуляція польоту в реальному часі з обраним азимутом і відхиленням вектора тяги, що вкладає ракету на необхідну траєкторію, яку ми перед цим перевірили. Логи польоту на даний час виводяться в консоль всередині середовища розробки (Qt Creator) та в CSV-файл для перегляду в Excel. Для ілюстрації деякі дані польоту екстраговано в таблицю з кроком 10 секунд, де зазначені основні події.

Приклад симуляції польоту

Початкові дані

Відстань до цілі: 662428 м

Дальність польоту з матриці чутливості по азимуту цілі 665458 м, відхилення вектора тяги двигуна 8,2°.

З корекцією азимуту, яку ми застосували, прогнозована дальність польоту 665315 м. Різницю опрацює система автоматичного маневрування. Це запас енергії для маневру.

Точка пуску: 50° 29′ 38,0″ N, 32° 58′ 35,2″ E

Точка цілі: 55° 45′ 13,5″ N, 37° 37′ 11,5″ E

Азимут на ціль: 0,456167 рад, або 26,1365°

Таблиця телеметрії

Похибка влучання: 6,6 м (X: 4,512; Y: —4,776)

Різниця з координатами цілі 0,2 кутові секунди.

При польоті на 662 км система автоматично скоригувала вектор тяги на старті, врахувала зміну щільності атмосфери з висотою від 0 до 151 км, прорахувала корекцію, зманеврувала на ціль і забезпечила влучання з похибкою лише 6,6 метри.

Висновок

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

Інженерам. Якщо ввести власні аеродинамічні залежності, масу, характеристики двигуна, то ним можна тестувати алгоритми керування під час пуску та наведення в процесі термінальної корекції, отримати дані по самим екстремальним режимам польоту.

Цивільним. Чисельно побачити як літає балістична ракета, дізнатись чому вона така небезпечна. Її швидкість, час прильоту і траєкторія, коли вона, як метеорит, просто падає з неба, дуже ускладнює перехоплення силами протиповітряної оборони.

Підписуйтеся на WhatsApp-канал DefTech спільноти!