МАТЕМАТИЧЕСКОЕ И ПРОГРАММНОЕ ОБЕСПЕЧЕНИЕ АНАЛИЗА РАССЕИВАНИЯ ТОЧЕК ПАДЕНИЯ ФРАГМЕНТОВ ЛЕТАТЕЛЬНЫХ АППАРАТОВ Текст научной статьи по специальности «Физика»

Электронный журнал «Труды МАИ». Выпуск № 61

www.mai.ru/science/trudy/

УДК 521.172:629.78

Математическое и программное обеспечение анализа рассеивания точек падения фрагментов летательных аппаратов

В.И. Гончаренко, Ю.С. Кан, А.А. Травин

Аннотация

Разработано математическое и программное обеспечение (МПО) анализа рассеивания фрагментов летательных аппаратов (ЛА) на поверхности Земли с учетом атмосферного торможения и малых случайных возмущений вектора начальной скорости на пассивном участке траектории (ПУТ). В качестве характеристики рассеивания выступает круговое вероятное отклонение (КВО) точки падения фрагмента ЛА от предполагаемой номинальной точки, определяемой без учета случайных возмущений. На модельном примере с использованием данного МПО показано, что КВО резко нарастает для номинальных траекторий с малыми начальными углами траектории в широком диапазоне значений сферической дальности полета.

Ключевые слова

летательный аппарат; возмущенное движение; торможение в атмосфере; баллистические производные; круговое вероятное отклонение

1. Введение

Актуальность исследования обусловлена в первую очередь тем, что решения по организации противодействия возможным угрозам для любой высокоразвитой страны принимаются на основе непрерывного мониторинга политической и ракетно-космической обстановки в мире или отдельных регионах Земли. С этой целью необходимо оценивать возможности множества ракетно-космических средств и систем, дислоцирующихся в различных позиционных районах Земли. Такие оценки базируются на решении баллистических оптимизационных задач определения зон досягаемости ЛА, зачастую в условиях неполного знания ряда проектных характеристик оцениваемых средств [1]. В таких условиях разрабатывать специальное математическое и программное обеспечение (СМПО) для решения задач оценки досягаемости традиционными способами, т.е. путём привлечения

большого числа разных специалистов, представляется весьма нерациональным. В то же время от степени точности оценок сложившихся в мире стратегических ситуаций существенно зависят расходы страны на создание новых и модернизацию существующих ЛА и средств противоракетной обороны [2].

Учет свойств области досягаемости имеет достаточно глубоко проработанную базу методического обоснования решений практических задач, возникающих при оценке угроз, связанных с применением различных ЛА [3,4]. Вопрос же учета влияния рассеивания отделяемых фрагментов ЛА исследован в меньшей степени.

При обосновании требований к параметрам объектов наземной инфраструктуры информационных систем слежения за полетом ЛА и обосновании требований к ней одним из важнейших параметров, используемых для учета влияния рассеивания на результаты пуска является КВО, характеризующее степень рассеивания точек падения на поверхность Земли. КВО является мерой кучности пусков при круговом рассеивании ЛА и представляет собой радиус круга, вероятность попадания в который равна 0,5 при условии совмещения центра нормального закона распределения ошибок пуска с центром круга. При нормальном законе рассеивания фрагмента связь между КВО, средним квадратическим отклонением (СКО) и срединным отклонением (СО) следующая КВО=1,1774 СКО=1,7456 СО [5]. Таким образом, в соответствии с терминологией [6], КВО представляет собой квантильный критерий качества надежности 0,5 для функции потерь, равной величине случайного отклонения точки падения фрагмента от центра нормального закона распределения. При этом КВО обычно бывает известным для типовых траекторий, например для траекторий максимальной дальности. Практический интерес представляет задача анализа зависимости КВО от свойств траекторий подлета ЛА с целью обоснованного выбора вышеупомянутых параметров.

В настоящей статье предложен подход к разработке СМПО, позволяющего провести расчёты по получению зависимости КВО от угла наклона траектории в начале ПУТ и полной сферической дальности полета ЛА, включающего помимо ПУТ также и АУТ (активный участок траектории). В качестве отправной информации принимается заданное соответствующее КВО, известное для оптимальных траекторий максимальной дальности.

Работа выполнена при поддержке Российского фонда фундаментальных исследований (грант № 12-08-00453-а) и ФЦП «Научные и научно-педагогические кадры инновационной России на 2009-2013 годы» (гос. контракт 16.740.11.0555).

2. Постановка задачи

Основные допущения следующие. Земля предполагается сферической. Вращение Земли не учитывается. Влияние атмосферы на АУТ не учитывается. Влияние атмосферы на ПУТ учитывается в соответствии с моделью торможения для фрагмента конусовидной формы со скругленной вершиной. Граница плотных слоев атмосферы принимается равной 90 км. Случайные возмущения, приводящие к рассеиванию точек падения фрагментов ЛА, моделируются нормальным распределением разброса вектора скорости ЛА в начале ПУТ с нулевым математическим ожиданием и ковариационной матрицей а2К0, где матрица К0

задана, а а2 - скалярный параметр. Распределение указанных случайных возмущений полагается одинаковым для всех допустимых траекторий ЛА. Величина ктах КВО для траектории максимальной дальности известна.

По заданным значениям параметров К0 и ктах требуется определить зависимость

КВО точки падения фрагмента ЛА от полной сферической дальности полета и угла наклона траектории ЛА в начале ПУТ.

В качестве исходных данных для модельных расчетов использовались характеристики ЛА Трайдент II, заимствованные из [7]. Матрица К0 в расчетах принята единичной, что

соответствует сферической модели рассеивания по скорости отделения фрагмента в конце АУТ. Но описанное ниже СМПО учитывает произвольный характер этой матрицы.

Величина ктах принята равной 100 ед. Это приводит к определению искомой

зависимости КВО в процентах от ктах. Например, при КВО=145 истинное значение искомого КВО рассчитывается по формуле КВО = 1.45 ХК^ ,

где ктах выражено в метрах.

3. Алгоритм оценки КВО

Для удобства обозначений принято: нижний индекс «0» у некоторого кинематического параметра означает, что этот параметр берется на начало ПУТ, а нижний индекс «а» говорит о том, что это параметр относится к моменту входа в плотные слои атмосферы.

В расчетах используются следующие вспомогательные системы координат (СК): х0, у, ^ - произвольная инерциальная, в которой задаются данные на начало ПУТ. В расчетах использована абсолютная геоцентрическая СК (АГСК).

ха, уа, ^ - инерциальная, связана с точкой Оа пересечения невозмущенной траектории фрагмента ЛА с границей плотных слоев атмосферы. Ось Оауа направлена по внешней нормали к границе атмосферы, ось Оаха ортогональна оси Оауа и направлена в плоскости невозмущенной орбиты в направлении движения фрагмента ЛА. Ось Оага дополняет систему до правой. Орты осей этой СК определяются по формулам

е1 - е2 Х е3

е2 -

(1)

ез -

V Х г

а а

V Х г

а а

где г - радиус-вектор точки Оа, ^ - вектор скорости фрагмента ЛА в точке Оа на невозмущенной траектории.

В, £ - плоская СК на поверхности Земли, связанная с точкой падения фрагмента ЛА при невозмущенном движении. £ - отклонение по дальности в плоскости невозмущенной орбиты, В - отклонение по боку.

Основные обозначения.

V - (Vх, vy, V ) - вектор скорости.

г — (х, у, 2 )Г - вектор положения.

ва — агоЩ (ууа /) - угол входа в атмосферу.

Ь — Ь (| V |, —ва ) - путь по Земле при торможении фрагмента в атмосфере.

4 — дЬ / д\>а\, Ьв—дЬ / д(—0а ) .

Эти производные определяются численно по трем точкам с помощью метода наименьших квадратов (МНК) для уменьшения влияния ошибок, обусловленных погрешностью численного интегрирования уравнений движения фрагмента ЛА в атмосфере. Процедура МНК для оценивания производных описана ниже в разделе 4.

а — (х , х , V , V , V ) .

"а у а? а? ха? уа? ха )

Я — (£, В )т.

г

а

С = дд / ду0 - матрица баллистических производных размера 5 х 3 . Вычисление этой

матрицы описано ниже в разделе 5.

А = дQ / дд - матрица частных производных размера 2 х 5 . Элементы матрицы А . Первая строка:

д, /а* = i, д, /= о, д, = /8v = l coso + ьа tg(0a)

11 а '12 а '13 xa v a 0 9

(1 + tg2 (0a ))| va\cos 0а

А. = 8S/8vva = L sin0 -Lel-„, ' -, As = 8S/8v = 0.

14 ya v а О /л , 2 / ^ \\ /г 5 15 za

(1 + tg (0a )) Ы COS 0a

Вторая строка:

A21 =8B / 8xa = 0, A22 =8B / 8za = 1, A23 =8B / 8vxa = 0, A24 =8B / 8vya = 0,

Д5 =8B / 8vza = L

v cos0

^mn - наименьшее и наибольшее собственные значения (2 х 2) -матрицы

ACK0CTAT.

Разработанный алгоритм включает в себя следующие три этапа.

Этап 1. Подготовка исходных данных для вероятностного анализа в результате моделирования дискретных пучков траекторий ЛА на ПУТ, соответствующих различным задаваемым с шагом 500 км значениям сферической дальности полета и различным значениям угла бросания в начале ПУТ. Подробно этот этап описан в разделе 6. Моделирование проводится с помощью разработанной для этой цели программы PrTrass2. На этом же этапе используется программа PassAtmProject для определения значений функции L(|vj,-0а) для смоделированного набора значений скорости и угла входа в атмосферу.

Этап 2. Определение параметров ,Лтх для смоделированных на этапе 1

траекторий. Этот этап подробно описан в разделе 5 и реализован в виде программы Ballist.

Этап 3. Для каждой смоделированной траектории искомое КВО к определяется по формуле:

ic = al

^"□f (у), (2)

где f (у) - КВО для двумерного нормального закона с плотностью

-11 ■

Г Л

p (xy )=xp 1" 2

2 y

x + —

л

а у = . Для нахождения f (у) в [6] разработана специальная процедура ее вычисления с

Лтах

заданной точностью, реализованная как функция программы Probably.

Параметр а определяется по известному КВО ^тах для траектории максимальной дальности. Для этого сначала выполняются первые два этапа алгоритма (программа PrTrass2 позволяет получать оптимальные траектории). Затем вычисляем

а = тах-. (3)

у Лтах f (у)

4. Процедура МНК оценки производной функции

Для численной оценки ряда частных производных используется следующий численный метод, основанный на МНК. Пусть требуется оценить производную f' (x0)

функции f (x) в точке x0, имея возможность лишь вычислять значения f (x) в различных точках xk с некоторыми ошибками (не обязательно случайными). Вычисляем значения функции в трех точках

У1 = f (x0 " h)+&, Уо = f (x0 ) + 4 , y2 = f (x0 + h) + 4

и с помощью МНК проводим на плоскости (x, y) прямую y = ax + b через точки (, yk ),

k = 0,2. Получившееся при этом значение параметра a и является оценкой искомой производной:

2 Y 2 2

Е ^ -1Е x хЕ yk

_ k=0_3 к=0_к=0

a — 0 .

2 2 1 Г 2 12 Е x2 Е xk

k=0 3 V k=0

5. Оценивание матрицы баллистических производных

В этом разделе описывается алгоритм определения матрицы C баллистических производных. Для реализации алгоритма раздела 3 необходимо рассчитать эту матрицу в СК

ха, у , ^. Обозначим и - вектор положения и скорости фрагмента при входе в атмосферу. Сначала определим матрицы

т =

r dvn

dxvh dXvh dXvh ^ (dvxvh dVxvh dVxvh 1

dvx 0 dVy 0 dVz0 dvx 0 dVy 0 dVz 0

fyvh fyvh fyvh J _ dVvh _ dVyvh dVyvh dVyvh

dVx 0 dVy 0 dVz0 ' V dv0 dVx 0 dVy 0 dVz 0

dZvh dZvh dZvh dVzvh dVzvh dZzvh

dVx 0 dVy 0 dVz0 J 0 dVy 0 dVz 0 J

в АГСК х0, у, ^. Это описано в подразделах 5.1 - 5.3. Затем вычисляем эти матрицы в СК ха, у , ^ в соответствии с правилом преобразования тензора 2-го ранга:

T = ТТ [ТТ T = ТТ [ТТ

га г уа V '

где строками матрицы Т являются вектора ег, е2, е3 (см. раздел 3). Искомая матрица С является составной:

( т Л

C =

T

1va\

T

У va3 J

где Туа1 и Туа3 - 1-я и 3-я строки матрицы Туа ■

5.1. Дифференцирование в АГСК по у0

Процедура определения производных основана на дифференцировании следующих

инвариантов задачи: векторного интеграла площадей

r0 Х v0 = rvh Х Vvh

(три уравнения), условия входа в атмосферу

. 2

= const = (R + 90km) (одно уравнение), где R - радиус сферической Земли, и интеграла энергии

(4)

(5)

vh

-П

(I 3 )=^-П(1 )

2 V у" ' 2 ,|и" (6)

(одно уравнение), где П - силовая функция центрального поля тяготения Земли.

Дифференцирование по у0 векторного соотношения (4) приводит к следующим трем соотношениям:

Уук + --^~Уу*к"2ук^— - 0 > (7)

^Ук д2ук

°>х0 ^хО

^хук дхук

°>х0 ^хО

^уук дуук

•л ¿Ук ^ Ук ^ уУ к Ук ^

°Ухо ^хО ^хО ^хО

дг^ Оу„,, дх„, ду.

ук , „ хук____ V ¿УЬ _ т ^Я^

У хук + 2ук ~ ~ У гук ХуЬ ~ - ¿0, (8)

дух0 дух0 дух0 дух0

дХк Л. дуук дУк _ дух*-- (а\

ууук + хук ~ ~ У хук у ук ^ - у0 • (9)

^хО ^хО ^хО дУх0

Дифференцируя по ух0 соотношение (5), получаем

Ох^+дгу*+Ог** - о • (1о)

дух0 ^хО ^хО

Дифференцируя по ух0 соотношение (6), с учетом (5) имеем:

д^хук , дууь , дугук

дух0 ^хО ^хО

Таким образом, получены 5 линейных уравнений (7) - (11) относительно шести неизвестных

дхУк ду ук дг ^ дУу дууук ду^^ Ох*

производных ——, ——, ——, ——, ——, ——. Первую из них, т.е. ——, определяем

дух0 д у х 0 дух0 ^хО д у х 0 ^хО д у х О

численно по методике раздела 4. Обозначим ее значение через а . Тогда для остальных пяти производных получаем 5 линейных уравнений:

Т~У2ук + уук"Г--Ч~ууук~2ук^-- 0 > (12)

дг, ду , ду ,

ук I _ хук „ гук _ _ , /1

Т—ухук + ¿ук --хук -Г— - ¿0 + ауук > (13)

^хО дУх0 ^хО

^уук дуук ду к пл,

- "у " ау уук > (14)

^хО дУх0 ^хО

, дг ук

у Ук Ук - "ахк, (15)

^хО д Ух О

Ухук + + ТТУкУ2Ук - Ух0 • (16)

х к у к

^хО ^хО ^хО

Правые части уравнений (12) - (16) обозначим посредством . Полученная система пяти линейных уравнений решается аналитически следующим образом. Из (13) получаем

дгА 1

"ук _

^хО Ухук

( ^хук , „ дУукЛ

Ь " 2 _— + X

Ь2 2ук ~ ^ хук ~ V ^хО ^хО J

А из (14) имеем

ук

8у,

х0

хук

Ъ3 - Хук

ЭУУУ^л, ^хук + Уук '

8у.

х0

8у.

(18)

х0 У

Подставляя (17) и (18) в (12), (15) и (16), получаем систему из трех линейных уравнений

V у хук

-Уук

хук Уу2к

8у.

хук

8ух у 2 ^

уук

+

уук

V ухук

хук 'ук

л8у . с

уук

У

8у„

Уук

V

Ук , хук ухук У

8ухук , Уукхук 8ууук , хук'ук 8у'

уу

V хук хук у

8у.

х 0

+ -

ухук 8ух0 ухук

уук

8у

ду

хук

8у

хук

8у„

+ у

уук

уук

8у„

+ у

х0

ду,

уук _ у уук

У 8ух0 у хук

у хук уук ^ у хук

ъъ + Уу^ъ2,

ухук

уук

уук

ду

= ъ

(19)

относительно трех неизвестных производных

^ук Эу

уук

8

уук

8ух0 Эух0 8ух0

. Система трех линейных

уравнений решается аналитически, поскольку обратная матрица к матрице размера 3 х 3

Эу и Эу ,

аналитически выписывается в соответствием с правилом Крамера. Производные ——, ——

Эух0 Эух0

определяются после этого из (17) и (18).

5.2. Дифференцирование в АГСК по уу0

8х. Эу .

Соотношения для определения производных

8хук дуук Эуук &хук Эууук Эуук

^у 0 ' ^у 0 ' ^0 ' Эуу 0 ' Эуу 0 ' ^у 0

выводятся аналогично подразделу 5.1. Производную

ду

ук

Эу.

определяем численно по методике

у 0

раздела 4. Обозначим ее значение через а. Тогда для производных

8хук

ук

8уу0 ^

получаются

у 0

следующие соотношения:

ук

1

С

8у.

у0

V,

Эх,

уук (

Ъ1 + у ук

^уук дУуук

уук

\

8у.

у0

8у.

ук

8уу0 ууук

Ъз+уук^

8у.

у0 \

V

8у,

- х

уук

к

у0

8у,

у0 у

где Ъ1 = - у0 - ауук , Ъ3 = х0 + аухук •

Система же трех линейных уравнений относительно производных

8у„,и Эу

"у 0

ду„Г1 ду.

у 0

(20) (21)

хук ук

ду

уук

ду.

у 0

получается следующей:

У

У

уук

4

уук у ук

уук

8у .

хук

ду

+

у0

^ук 8у:

У ^к ^ ухук _ -хук--у

хук

Ууук Эуу0

у

V уук

ук

у

V уук

ук

уук

8у

уук

у

ду.

+

у0

уук

у

V уук

уук- х

ук

8у

уук

у2 ]дууук_ + уу_куук 8у-.

"ук

уук

уук у

ду.

у 0

ууук 8уу0

= Ъ4 +

ду.

уу

+ -укЪ3,

у0

уук

ук

уук

Ъ -—Ъ

ууук

31

ду.

хук

ду.

хук

ду.

+ у

уук

уук

у0

ду

+ у

ду.

уук

уук

у0

ду.

= Ъ

у0

уук

где Ъ4 =-ауук , Ъ5 = уу 0 •

5.3. Дифференцирование в АГСК по у0

Эхук ^ук Эуук духук 8ууук Эу

Соотношения для определения производных

уук

Эуу0 Эуу0 Эу у 0 8у у 0 Эуу 0 Эуу 0

Эу ь

выводятся аналогично подразделу 5.1. Производную —— определяем численно по методике

ду„п

раздела 4. Обозначим ее значение через а . Тогда для производных

Эуук дхук

Эуу0 8уу0

получаются

следующие соотношения:

су^

(

ук

8уу0 уук

Ъ1 - уук

8у.

у к

ду„,

+ у.

к

дх.

ук

8уу0 ууук

8ууук 1 8у у0 у

8у

г. _ ^ хук , ■.„ уук

Ъ2 уук ~ + хук

ду

8у,

у 0

8у.

у 0 у

(23)

(24)

где Ъ1 = у0 + аууук , Ъ2 = -х0 - ауук •

- 8ухук ^уук ду ук

Система же трех линейных уравнений относительно производных ——, ——, ——

8уу0 8уу0 Э у у 0

получается следующей:

у к

V ,

V уук

л

уук у ук

ду.

(

х к

у

ду

+

у0

хук

V

ххк

уук ^уук у

ду.

у к

ду.

+

у0

хук

v ,

V уук

уук

Ук у хук у уук у

8

у к

ду

V , V

у0

у к

у к

Хук^ук ду хук | Уукук 8ук

гх2 х,

ууук 8уу0

V.

у к

ду.

у0

ук уук

v ,

V уук

уук у

^ уук _ Ъ | ук Ъ уук А

Ъ4 + Ъ2 Ъ1,

ду

у0

у к

у к

ду

х к

ду

х к

ду

+ у

у к

у к

Эу„

+ у

ду.

у к

у к

ду,

=Ъ

(25)

1

где Ъ4 =-ауук > Ъ5 = уу0

6. Моделирование траекторий

6.1. Расчет начальных условий баллистического полета на заданную дальность

Отметим, что поставлена задача исследования некоторых свойств баллистического полета (БП) с использованием сферической модели Земли без учета влияния атмосферы и вращения Земли. Исходными данными являются требуемая сферическая дальность полета /сф и определенные ниже в подразделе 6.3 параметры конца АУТ - дальность и высота полета на АУТ: /к и кк. Требуется определить в векторной форме радиус вектор и скорость полета на момент начала БП, а также координаты точек падения в целевой системе координат, соответствующие выбранным начальным условиям (НУ). Предполагается, что величина и наклон вектора скорости на момент окончания АУТ будут варьироваться с помощью датчика случайных чисел в определенных пределах.

В сферической системе координат дальность баллистического полета /б определяется параметрами конца АУТ, выступающими в качестве НУ БП: /б = /б(Кк, вк, Лк), где 0к - наклон вектора скорости Ук в конце АУТ (в точке К); кк - высота конечной точки АУТ. Если требуется определить значения перечисленных параметров для случая полета на заданную дальность (что требуется в настоящем исследовании), то возникает краевая задача, имеющая аналитическое решение (исчерпывающе подробно эта задача исследована Д.А. Погореловым [9]). С учетом того обстоятельства, что в данной задаче /к, - известная дальность полета на АУТ, требуемая дальность БП /б и соответствующая ей сферическая дальность (Р) БП вычисляются по формулам:

/б = /сф - /к, Р = /б / Я (26)

Здесь радиус Земли Яз = 6371 км.

Уравнение, связывающее НУ БП при заданной угловой дальности Р (уравнение годографа скорости) имеет вид [9]:

V2 - —-2-—, (27)

sin 20к ctg F + cos2 9 ^2 + h + h ctg2 F

где относительная высота полета (параметр h) рассчитывается по формуле: h

R + h

3 2

гравитационный параметр b0 = 398600.44 км /с , а V - требуемая для полета на заданную угловую дальность скорость в начале БП при заданном же угле наклона этой скорости.

Основным вариантом расчетов при использовании рассматриваемого алгоритма является случай, когда параметр вк задается в исходных данных (варьируется от минимального до максимального). В этом случае алгоритм задачи расчета НУ исчерпывается формулами (26) - (28).

Второй вариант решения задачи расчета рассматриваемых начальных условий заключается в получении их из условия расчета минимальной начальной скорости полета на заданную угловую дальность Р. Согласно [9] алгоритм описывается формулами (29), (30) с учетом (28):

( ^ Г Л

„ 1

в = — агс сог к 2

F к +-

^ 2 ^п F

(29)

F

2 (30)

V - „

г

Таким образом, формулы (29), (30) замыкают алгоритм расчета НУ БП.

Для исследования влияния вариаций НУ БП на рассеивание точек падения требуется пересчитать параметры скорости и положения центра масс ЛА в точке окончания полета на АУТ из плоской сферической системы координат (СК) в пространственную СК. Для удобства решения рассматриваемой задачи на этапе расчета координат точек падения целесообразно использовать «абсолютную геоцентрическую систему координат» (при учете вращения Земли - «относительную геоцентрическую систему координат»).

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

Алгоритм пересчета НУ БП из используемой выше сферической СК в АГСК приводится ниже.

Хн0 -

- соб А СОБ 4 ^п В - А ^п 4

Бт А СОВ 4 - СОВ А ^П В 4 СОБ А СОБ В0

(31)

>0 -

СОБВ СОБ 4 соб В0 Б1П 4 Бт Вп

(32)

1но - Хно X ^ . (33)

Формулы (31) - (33) служат для пересчета параметров, заданных в проекциях на оси НГСК в АГСК. В них: Хн0 - орт оси Охн спроектированный на оси АГСК, ун0 - орт оси Оун

спроектированный на оси АГСК, 1нй - орт оси Охн спроектированный на оси АГСК.

Непосредственное вычисление НУ БП с учетом приведенных формул выполняется тривиально:

Г0 - Хн0 зт 4 + ун0 соб(34)

(35)

I

Т° - Хн0 С0Б4 - Ун0 ^П 4

- Г0 (Я + кк), а 4 -

Гк - г

Яз

V -(

-(Г0 ътв„+т° соб

(36)

(37)

Возмущенный вектор скорости в момент начала полета на ПУТ вычисляется по формулам:

Vвозм - V + АV

(38)

Здесь А¥ -

А^,

, где А¥х, АУу - случайные вариации соответствующих компонент

вектора А V , (АУ2 = 0, т.е. боковые возмущения не учитываются) подчиненные нормальному закону распределения и имеющие нулевое математическое ожидание и задаваемое среднеквадратическое отклонение, а лх, ку - признаки формирования вариаций соответствующих компонент скорости (признаки задаются единицей или нулем, соответственно).

6.2. Расчет параметров в точке падения, соответствующей выбранным НУ

Расчет координат точки падения и скорости полета в ней, соответствующих НУ V и Г (как в случае номинального, так и для возмущенного полета на ПУТ), осуществляется по алгоритму, представленному ниже.

V2 = (/7возм )т Кквозм; Уквозм г к - Яз + Ик, У

2Ъ

(39)

лвозм и

= , ¡1 -

( Vв°змр V

у к к

V возм Г

V к 'к у

е2 = 1 - 4 (1 -у) 6»во

_ _ тт-возм „„ лв

а - V Гк СОБ^.

V = V 2 - 2Ъ0(к- к)

1'' гк (л+А„) '

0„ = - arccos

а

V (л + к)

^возм = 2 arctg

у cos2 0ВО

- 2 - К, К =

г - л - К

Л + Ац

Гц° = :е„0 sin ^возм + у«0 C0s ^возм , Гц0 = хрНо COs ^возм - у«0 sin ^возм ,

г = г

ц ц

Л, V =( г0 ^п6>ц + ?0^6>ц )КП

г V2

ц ' 2

2Ъ

(40)

(41)

(42)

(43)

(44)

(45)

(46)

(47)

Г и V - координаты и скорость полета в точке падения в АГСК.

Удобный для реализации на ЭВМ алгоритм расчета времени полета на пассивном участке траектории (имеется и альтернативный вариант) также описан в [9]:

~ 1 - 2у

Д =-'-; у = arcsin Д; Д = е ^ у; г = 1,2.

е

/1 = 1 -2У1, / = 1 -2У2, у =7е2 -./12, У2 Р = -, Рг = -

* * е "

Л

е

т =

У

а

ж £

—+ у- аг^ — + ц

2 У2

ж А

—+ у - аг^ — 2 У

Л"

у

, где ц = sign ^т0возм ).

(48)

(49)

(50)

1

Ъ

0

6.3. Проектировочный баллистический расчет

Полная дальность полета складывается из дальности активного (/к) и пассивного (Ьп) участков траектории (рис. 1) [9, с.86-89].

Рис. 1. Траектория движения и схема сил, действующих на ЛА

Полет на ПУТ состоит из движения в безвоздушном пространстве (/эл) и на атмосферном участке (/а). В дальнейшем изменением дальности за счет сопротивления воздуха на атмосферном участке пренебрегаем, и весь пассивный участок считается по эллиптической теории.

Известно, что дальность Ьа определяется параметрами конца АУТ: Ьп = Ьп(¥к, 9к, Лк), где 9к - наклон вектора скорости Ук в конце активного участка (точка К); кк - его высота. Тогда полная дальность полета: Ь = Ь( ¥к, 9к, Ик, /к). Эта же зависимость при заданном угле 9к (т.е. при выбранной программе движения на АУТ) может быть записана через основные проектные параметры (ОПП), которые в процессе проектирования должны быть выбраны так, чтобы выполнялось ограничение Ь > Ьтах, являющееся условием обеспечения заданной дальности полета.

Управляемые ЛА одного класса имеют обычно подобные одна другой программы движения на АУТ. Это приводит к тому, что величины 9к, Ик, /к для различных ЛА с одинаковой полной дальностью полета оказываются примерно постоянными. Величины /к и кк пропорциональны коэффициенту тяговооруженности, поэтому их значения следует умножить на поправочный коэффициент К = ^ср/0.5, где ^ср - среднеинтегральная величина коэффициента начальной тяговооруженности для всех ступеней.

При известных величинах Ик, 1к и вк скорость Ук ЛА в конце АУТ однозначно зависит от дальности полета. Для расчета величины Ук по эллиптической теории используются следующие рабочие формулы:

Ё = ^, (51)

а = ага%-,-^^ПЁ-, (52)

(Я + К)" Я соё

V =

1

Ь0 ^ (рс а

^ —„ ; „ч, (53)

Я + К 2 ео8 в соб (а-вХ

з к к \ к /

где Кз = 6371 км - радиус Земли; ¿0 (параметр гравитационного поля Земли (¿о = 398600.44

3 2

км /с , - постоянная величина). Скорость в конце АУТ:

п 1

V =Е-^поТ, (54)

;=1 1 Му

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

Гк+АГпот= ку-Гк, (55)

где ку - коэффициент потерь скорости, зависящий от дальности полета, удельной тяги и начальной тяговооруженности ступеней. Для ЛА с дальностью полета Ь - 10... 14 тыс. км коэффициент ку = 1,15... 1,25, при этом большим дальностям полета соответствует меньшее значение kу.

Если приближенно принять, что удельные тяги двигателей всех ступеней одинаковы, то формула (55) примет вид:

¥К = Лд.ср. -„ V л-- - А^от, (56)

у р (1 -Мт1 )(1 -Мт2 )-"(1 -Мтп )

где ./уд ср - среднее значение удельной тяги двигателей всех ступеней. Для ЛА с ЖРД можно приближенно принять:

1 (1 +1 п 1 1 1 уд.01 +1 уд.п1 т

1удср- 2п -1

+ 1уд.п1 (57)

г=2 у

Для ЛА с РДТТ:

1 ( 1 + 1 ' =п 1

1 1 уд.01 + 1 уд.п1 ж—I т

1 уд.ср

п

+ Е 1уд. п, (58)

г=2 у •

Частная производная Ьу используется для коррекции дальности полета при уточнении

коэффициентов заполнения топливом ступеней (ц-ц). Сначала вычисляется Д«пр - ———,

^У ' Уцпр

дV

где V = ——, а ЛЬ - расхождение между требуемой дальностью полета и реализуемой при

ЭЦпр

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

7. Расчет торможения в атмосфере

Процесс атмосферного торможения существенно зависит от аэродинамических свойств отделяемого фрагмента ЛА. В настоящей статье в качестве фрагмента рассматривается осесимметричное твердое тело конической формы со скругленной

Рис. 2. Геометрические параметры фрагмента ЛА Интегрируются уравнения движения в геоцентрической, вращающейся вместе с Землёй, системе координат

г = V V = 2

вр вр ей

где г - радиус-вектор центра масс фрагмента. Полное ускорение 2П фрагмента в предположении нулевого значения угла атаки может быть оценено равенством

2 = 2(ет) - 2( ) - Сх х(а) п т вр

2(ю) - кориолисово ускорение,

2(ю) =-2[й3 х увр ] - [Й3 х [Й3 х г]] - угловая скорость вращения Земли (в расчетах в силу принятых допущений принята равной нулю), 2( гг) - гравитационное ускорение. Для Ньютонова поля

g gr) = g( gr) = JL f R = 1

g gH R3 '

8 - площадь миделя, т - масса фрагмента, двр - скоростной напор, Сх - коэффициент

лобового сопротивления.

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

Геометрию боевого блока опишем параметрами г, ^, ^ и Д, см. рис. 2. Ударную

составляющую коэффициента лобового сопротивления возьмём в соответствии с примером [10]:

Cx, a Ро

2 \ cos2

\cos2 Дк

R

-d

V

+ sin3 £(2R -hsin Д) Ra

м

где p0 - коэффициент давления в точке полного торможения, определяемый, например опытным путём.

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

8. Программное обеспечение

Разработанный программный комплекс включает в себя имеющие самостоятельное значение программы PrTrass2, Ballist, Probably и программу PassAtmProject расчета торможения отделяемого фрагмента в атмосфере. Для использования данного СМПО

специальные требования отсутствуют. Достаточно наличие персонального компьютера с операционной системой Windows любой версии. Для расчета искомой зависимости КВО от кинематических параметров начала АУТ указанные программы многократно запускаются в следующем порядке: PrTrass2, PassAtmProject, Ballist, Probably. Значения угла бросания и дальность полета перебираются по узлам плоской дискретной сетки возможных значений.

8.1. Описание и порядок выполнения программы PrTrass2

Программа предназначена для расчета условий входа в атмосферу при высоте Ha=90 км для случая траектории минимальной энергии, соответствующей задаваемой дальности.

Предусмотрена возможность моделировать полет ЛА с разной тяговооруженностью. При этом коэффициент тяговооруженности и требуемая сферическая дальность задаются в исходных данных (файл «TrassID.txt»). Их можно скорректировать и после запуска программы и это будет учтено в расчетах. В этом же файле можно задать значение «по умолчанию» для остальных параметров, задаваемых на главной форме программы. Изменение этих значений в процессе работы программы, т.е. после запуска программы щелчком (одинарным или двойным, - в зависимости от настроек операционной системы), -повлияет на их значения в текущем сеансе работы, но при повторном запуске все значения управляющих параметров вернутся к исходным значениям, заданным в файле «TrassID.txt». Изменение любого параметра непосредственно в указанном файле скорректирует значение этого параметра по умолчанию для всех последующих запусков программы.

Имеется возможность дополнительно выдать на печать «трассу полета», параметры: t, B, L, H с «шагом печати» по времени. На значение этого параметра по умолчанию распространяется все сказанное выше в отношении всех остальных параметров.

Основные результаты расчетов печатаются в файл «Rez.txt», который после нажатия кнопки «Выход» автоматически открывается в компоненте Windows, заданном в системных настройках для открытия текстовых файлов.

Для корректной работы программы в операционной системе должны быть заданы русские региональные параметры (т.е. разделителем целой и дробной части числа с плавающей запятой задается точка, а не запятая, которую OS Windows предлагает «по умолчанию»).

Для запуска программы без компиляции (рекомендуемый режим работы) необходимо два раза (или один, если так задано в настройках OS Windows) по выполняемому файлу

«PrTrass2.exe», который отличается от всех остальных «иконкой», на которой изображен земной шар.

Все выходные параметры программы имеют интуитивно понятные (типичные для этих параметров) обозначения и не требуют дополнительных пояснений.

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

траектории для выбранной дальности (т.е. автоматического расчета в°р) необходимо щелчком манипулятора «мышь» рядом с надписью на панели управления вычислительным процессом «оптим-я. траектория» изменить режим расчета.

Имеется возможность задавать возмущения, которые прибавляются к вектору скорости в конечной точке АУТ в АГСК (в этой СК осуществляется и моделирование полета на ПУТ). При этом, возмущение задается пользователем в процентах от номинала той составляющей вектора скорости, которую выбирает пользователь, а его знак вне зависимости от знака компоненты задается пользователем.

Поскольку возмущение задается в точке начала ПУТ (корректируется вектор скорости бросания), а точка входа в атмосферу находится вблизи конца ПУТ, в точке входа изменяются все 3 составляющие вектора входа в атмосферу.

8.2. Описание и порядок выполнения программы Ба1Н81

Программа Ballist.exe предназначена для расчета баллистических производных и параметров полета фрагмента ЛА при его входе в атмосферу при высоте На = 90 км. Расчет производится в три этапа: расчет матрицы баллистических производных, расчет частных производных пути фрагмента по скорости и углу входа в атмосферу и расчет собственных значений матрицы...

Интерфейс программы изображен на рис. 3.

Рис. 3. Интерфейс программы Ballist.exe Входные данные в программу вводить не нужно, вся необходимая информация

извлекается из файла данных data1.xls. Необходимо следить за тем, чтобы данный файл

находился вместе с программой в одной директории. Исходными данными для программы

содержаться в таблице вариации модуля скоростей Ук и углов Theta в момент входа

фрагмента в атмосферу и соответствующие им значения пройденного пути, таблице

покомпонентной вариации вектора Ук и таблица вариаций в зависимости от угла Theta.

При изменении файла данных настоятельно рекомендуется не изменять формат значений

исходных параметров, а также структуру и расположение таблиц, т.к. координаты таблиц

жестко прописаны в программе.

Расчет параметров полета фрагмента может производиться для оптимальной

траектории или для другого заданного угла входа. При выборе угла для оптимальной

траектории в табличке «Что будем вычислять?» табличка «Выберите значение угла входа в

атмосферу Theta» остается пассивной, а в случае выбора определенного угла данная

табличка активизируется и пользователь может выбрать значение угла входа из

предложенных программой (см. рис. 4).

Расчет параметров

ШШ

File

Что буцем вычислять?

Произвести расчет только для оптимальной траектории

(• Произвести расчет для заданного угла вксща

Рассчитать матрицу баллистических производных С

Рассчитать частные производные для пути L

Вычислить собственные значения

Собственные значения матрицы А Баллистические производные

Lambda_nnin = dL_d^vh =

Lambda_nnaK = dL_dTtvh =

Gamma =

sqrt(Lambda_rriaK] =

Собственные значения матрицы Грамма

Lambda_1 = Lambda_2 =

Garnma2 = sqrt(Lambda_2) =

Рис. 4. Выбор угла входа в атмосферу

После выбора угла входа в атмосферу предлагается рассчитать матрицу баллистических производных векторов яуь и ууь. Это действие осуществляется с помощью кнопки «Рассчитать матрицу баллистических производных С». Сформированная матрица отображается в правом нижнем углу интерфейса (см. рис. 5).

Выберите значение угла входа в атмосферу Thêta (в грацусак) Г 15 Г 25 С @ Г 45 Г 55

Рис. 5. Расчет матрицы баллистических производных

На втором этапе расчета производится вычисление частных производных dL/dVk и dL/dTheta. Данное действие позволяет осуществить кнопка «Рассчитать частные производные для пути L». После выполнения процедуры значения производных приводятся под надписью «Баллистические производные» (см. рис. 6).

Рис. 6. Расчет частных производных для пути L

Кнопка «Вычислить собственные значения» выводит на экран собственные значения матрицы А и матрицы Грамма скалярных произведений векторов с частными производными dRvh/dVk для отклонений по x и y. Кроме того, на экран выводятся отношения минимального к максимальному собственных значений (Gamma и Gamma2) и корни квадратные из наибольших собственных значений (см. рис. 7).

Расчет параметров

File

Что буцем вычислять?

с П роизвести расчет только для оптимальной траектории Г* Произвести расчет для заданного угла б-ода

Рассчитать матрицу баллистических производным С

Выберите значение угла входа в атмосферу Thêta (в градусах) Г 15 С 25 (* 35 Г 45 Г 55

Рассчитать частные производные для пути L

Вычислить собственные значения

2-я строка RTenzor_ -0.0002331 -0.0001207 -0.0007211

Xk, Yk, Zk = 1031757.00000 2834728.00000 5787204.00000

Собственные значения матрицы А

Lambda_min = 17321 Э.3816817384 Lambda_maw = 12787495.4542238358

Gamma = 0.0135458975 sqrt(Lambda_rnax) = 3575.8807735857

Собственные значения матрицы Грамма

LambdaJ =147397.4873585788 Lambda_2 = 10133030.8843582779

Gamma2 = 0.0145462387 sqrt(Lambda_2) = 3183.2422095026

Баллистические производные

dL_dVvh = 1.4888888988 dL dTtvh = -258403.8858041388

-1041.7330233071 -539.5306081571 -2786.2331306805

-0.0014078585 -0.0007210840 1041.7330567820

1.1885010470 1.1315123531 1.5806243635

1.8388124387 -0.1077315881 2.2185781842

0.0000000087 0.0000000122 -0.0082208480

I

Рис. 7. Вычисление собственных значений

Выход из программы осуществляется с помощью опции «Close» главного меню «File» или нажатием красного крестика в верхнем правом углу программы.

8.3. Описание и порядок выполнения программы Probably

Программа предназначена для расчета:

• нижней и верхней оценок вероятности попадания нормально распределенного двумерного случайного вектора в круговую область заданного радиуса с центром в математическом ожидании;

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

Рассматриваемый двумерный вектор имеет нормальное распределение с нулевым

вектором математического ожидания и диагональной ковариационной матрицей с

2 „

дисперсиями 1 и о . Интерфейс программы представлен на рис. 8.

Расчет вероятности н кващщщ попадания двумерного случайного вектора в эллиптическую область

Рис. 8. Интерфейс программы Probably.exe

Исходными данными программы являются значение среднеквадратического отклонение второй компоненты случайного вектора, радиус круга, доверительная вероятность попадания и точности оценивания границ вероятности и квантили.

Если необходимо рассчитать границы вероятности, то исходными данными для этой операции являются значение среднеквадратического отклонение, радиус круга точность оценки; если же мы оцениваем квантиль радиуса области, то необходимыми данными будут среднеквадратического отклонение, доверительная вероятность и точность оценки границ квантили. Все эти данные могут вводиться и корректироваться в соответствующих полях и имеют формат десятичного числа с разделителем запятая. Для нахождения КВО уровень доверительной вероятности необходимо положить равным 0,5.

Действие программы осуществляется после нажатия кнопок «Оценивание вероятности» или «Оценивание квантили», в зависимости от того, что нужно оценить.

Выходными данными являются, соответственно, нижние и верхние границы вероятности и квантили, а также разности этих границ (см. рис. 9).

Расчет вероятности и квантили попадания двумерного случайного вектора в эллиптическую область Г-П!^

Рис. 9. Выходные данные программы

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

Выход из программы осуществляется с помощью красного крестика в верхнем левом углу интерфейса.

8.4. Описание и порядок выполнения программы Ра88А1шРго]ее1

Программа PassAtmProject моделирует атмосферное торможение фрагмента в соответствии с методикой раздела 7. При запуске программы открывается интерфейс, содержащий опции «Редактирование», «Просмотр результатов», «Запись», «Выход».

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

Расчет кинематических параметров при при торможении фрагмента осуществляется задействованием опции «Запись», где можно указать путь и имя текстового файла, куда будут записаны результаты. В этом файле будет записана таблица зависимостей параметров движения от времени. Для реализации алгоритма раздела 3 в этой таблице требуется лишь

последнее значение в колонке «Путь по Земле (км)», равное полному пути по Земле при движении фрагмента в атмосфере.

9. Результаты расчетов

2000 3000 4000

—1-1-'-1-1-1-'-1-'-'-1-'-1-1-1-1-1-1- Угол бросания

20 30 40 50

Рис. 10. Номограммы изменения КВО в зависимости от угла бросания

Результаты расчета КВО в зависимости от угла бросания с шагом 1000 км по дальности полета представлены на рис. 10. Можно отметить два обстоятельства. Во-первых, видно, что для «настильных» траекторий [5], характеризующихся малыми углами бросания, рассеивание фрагментов возрастает в разы и может оказаться неприемлемым. Это возрастание носит «взрывной» характер, о чем свидетельствует также табл. 1, в которой представлены результаты для дальности 8000 км. Сравнивая эти результаты с результатами, полученными с использованием предложенных моделей, но без учета аэродинамического торможения, можно сделать вывод о том, что отмеченный «взрывной» эффект является следствием именно аэродинамического торможения.

Таблица 1. Расчеты для дальности 8000 км.

0к Ы УуЪ \/Апах У / (у) К

15 15.88 6867.0 11453 0.0017 0.676 58.1

25 25.63 6683.6 4967 0.008 0.681 25.4

35 35.39 6754.2 4386 0.01 0.682 22.4

45 45.18 7096.3 4558 0.01 0.682 23.3

55 55.01 7804.3 5001 0.011 0.683 25.6

Во-вторых, можно заметить существование критического угла бросания, выше котрого функция КВО практически не меняется во всем диапазоне допустимых дальностей. В рассмотренном примере это критическое значение равно примерно 25 градусам. При углах бросания ниже этого критического значения рассеивание фрагментов резко возрастает. Более точное определение критического значения угла бросания требует использования более точных моделей движения ЛА и более точных моделей внешних факторов (форма Земли, свойства атмосферы).

10. Заключение

Разработано СМПО анализа рассеивания фрагментов ЛА на поверхности Земли с учетом атмосферного торможения и малых случайных возмущений вектора начальной скорости на ПУТ. В качестве характеристики рассеивания использовано КВО точки падения фрагмента ЛА от предполагаемой номинальной точки, определяемой без учета случайных возмущений.

Предложенное математическое и программное обеспечение позволяет при сделанных допущениях на качественном уровне анализировать потенциальные возможности различных ЛА с отделяемыми фрагментами конусовидной формы. На модельном примере обнаружено, что рассеивание фрагментов резко нарастает при использовании «настильных» траекторий. Данный эффект является следствием аэродинамического торможения при движении фрагмента в атмосфере.

Библиографический список

[1] Проектирование и испытания баллистических ракет / под ред. В.И. Варфоломеева, М.И. Копытова. М.: Воениздат, 1970.

[2] Фатеев В.Ф., Суханов С.А. Концепция развития РКО России // ВКО, 2007, №2 (33).

[3] Василенко В.В., Гончаренко В.И. Разработка программно-информационного комплекса для визуализации областей земной поверхности, досягаемых баллистическими ракетами // Вестник компьютерных и информационных технологий, 2010. № 7.

[4] Василенко В.В., Гончаренко В.И., Шестаков П.С. Построение специализированных картографических проекций и зон досягаемости баллистическими объектами // Вестник МАИ. - 2011, т.18, № 6.

[5] Военный энциклопедический словарь Ракетных войск стратегического назначения / Гл. ред.: И.Д.Сергеев, В.Н.Яковлев, Н.Е.Соловцов и др. М.: Большая Российская энциклопедия, 1998.

[6] Кибзун А.И., Кан Ю.С. Задачи стохастического программирования с вероятностными критериями. М.: ФИЗМАТЛИТ, 2009.

[7] Гронлунд Л., Райт Д. Баллистические ракеты морского базирования с настильными траекториями: техническая оценка и возможности контроля // Наука и всеобщая безопасность. Технические предпосылки для политических инициатив по контролю над вооружениями и проблемам окружающей среды. 1992, Т. 3, № 2.

[8] Кан Ю.С., Травин А.А. О приближенном вычислении квантильного критерия // Автоматика и телемеханика, 2013, № 5 (в печати).

[9] Погорелов Д.А. Теория кеплеровых движений летательных аппаратов. М.: Физматгиз,

1961.

[10] Краснов Н.Ф. Аэродинамика. Ч. II: Методы аэродинамического расчета. - изд. 3-е, перераб. и доп. М.: Высш. шк., 1980.

Сведения об авторах:

Гончаренко Владимир Иванович, д.т.н., доц., директор Военного института Московского авиационного института (национального исследовательского университета) (МАИ). Волоколамское ш., 4, Москва, 125993; тел. (499)158-17-01, e-mail: vladimirgonch@mail.ru

Кан Юрий Сергеевич, д.ф.-м.н., проф., профессор кафедры теории вероятностей Московского авиационного института (национального исследовательского университета) (МАИ).

Волоколамское ш., 4, Москва, 125993; тел. (499)158-41-13; e-mail: yu_kan@mail.ru

Травин Андрей Александрович, аспирант Московского авиационного института (национального исследовательского университета) (МАИ).

Волоколамское ш., 4, Москва, 125993; тел. (499)158-5G-13; e-mail: dron-mail2GG1@mail.ru