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

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

www.mai.ru/science/trudy/

УДК: 629.764::519.622

Приближенные методы расчета внеатмосферного активного

участка траектории

Карагодин В. В.

Институт стратегической стабильности Госкорпорации Росатом, Луганская, ул., 9, Москва, 115304, Москва, Россия e-mail: mikhailov@iss-atomt. ru

Аннотация

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

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

уравнения движения, внеатмосферный активный участок траектории, метод решения, точность, быстродействие

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

упрощения решения, не выражающегося через элементарные функции, в работах [1,2] вычисление неберущихся интегралов удерживает только квадратичные члены. В итоге, полученные расчетные формулы имеют невысокую точность на протяженных траекториях. Так, в соответствии с приведенным в [2] примером расчета траектории выведения второй ступени ЛА с начальной тяговооруженностью п0=0,8 и удельной тягой Руд=445 секунд из начальной точки (скорость 3073 м/с; высота 40 км; угол наклона вектора скорости к местному горизонту 20о) в конечную перигейную точку (высота 200 км, скорость 7918 м/с) методические ошибки достигают 56 км по координатам (относительная погрешность ~3%) и 120 м/с (1,5%) по скорости, хотя конечная масса ступени определяется с приемлемой погрешностью 0,3%. У метода из работы [1] погрешности еще более высокие - 150 км (8%) и 205 м/с (2,5%) соответственно. Для сравнения, при таких же затратах процессорного времени погрешность численного решения уравнений движения стандартным методом Рунге-Кутты за один шаг интегрирования (к «392 с) составляет 13 км (0,7%) по местоположению и 82 м/с (1%) по скорости. То есть в рассматриваемом примере только численный метод удовлетворяет допустимой в проектных расчетах погрешности определения скорости и координат, равной 1% [1].

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

Общее решение линеаризованных уравнений движения в проекциях на оси начальной стартовой системы координат при начальных условиях ^ = 0, х = х0, у = у0, V = Ко, К = Ко приведено [2]:

Fn . 1 г —cosp . x = x0 cos(01 л—sin g)0t л--I-sin(0(t -z)dz;

¿¡л ¿¡л J iiJ

(o (o 0 m

V \ t p [ft

y = yQch(Q42t л--sha0yp2t л---¡= I —sinpsh(0V2(t - z)dz - (ch(0V2t -1).

(0V2 (0V21 m 2(0

где (0 = ylg3 / R, - круговая частота (частота Шулера);

g3 - ускорение гравитационного поля на поверхности Земли; R - радиус Земли;

— = urn - сила тяги ракетного двигателя; u - пустотная эффективная скорость истечения; m = m - mt - масса ЛА;

p = p0+pt - угол тангажа.

Неберущийся интеграл в выражении для горизонтальной координаты x записывается через интегральные синусы Si и косинусы Ci:

г P cos p . . . , и . т т 1 т 1 т . I-sm^0(t-т)ат=— (aA + a-^2-bj -bj2),

о m 2

где a = sin(^ - pT ); a2 = sin(^ + pT ); b = cos(^ - pT); b2 = cos(^ + pT); а = Юо(т-1); (t =(o+(T; t =

m

T-t

= j cos(^0 ( dt = J T

7

t cos(^ + ()7

7

t sin(^

7

• sin(^ + Ф)7

Ci(rn0 - ()T - Ci(rn0 - ()(T -1) при ;

i т -1 .

при ( = 0°;

Ci(®0 + ()T - Ci(®0 + ()(T -1) при ( ф -ф0 ; Т -1

при ф = -Юо';

12 = J COS(^0 + ()7 dr =

т 7

= - sin(0° - ф)т dT = Si(p - ()(Т -1) - Si(p - ()T ;

T 7

J = J sin(0° +()r dT = Si(0o + ()(T -1) - Si(rn0 + ()T .

J T

T 7

Неберущийся интеграл в выражении для _у можно представить в виде: г Р sin( - T)d т = - U (a3/3 - a J4 - b3J3 - ) ;

J m 9

0

-ctt/2 • CTV- ■

где a = е sin(; a4 = е sin(;

г. -o4— т. aV-

b3 = е cos (pT; b4 = e cos (T; r г ег cosаг ^^ x Г ег sinac ^

I3 = J -7-dC; J3 = J ---dC ;

®°V-T c ®°V-t c

r г ег cosa^ ^^ x Г ег sinаг ^ (

h = J -г-d4 = J ""; a = ^7-.

-®0V-t г -®0 V—т c ^ov-

Разложив в подынтегральных выражениях тригонометрические функции при а < 1 (скорость разворота по тангажу ( < 0,1 град/с) или экспоненту ег при а > 1 в ряд Тейлора,

получим рекурсивный алгоритм вычисления интегралов I , I , J , J :

г £ ^ 2n г ' £

при a < 1 J= щ + Z(-1)п Je-SlnaCdг = ]Г (-1)й+^-0-^;Е2й-2;

c n=1 (-n)! c n=1 (-n 1)!

<

T

£ ш гп

Щ=\ -¿£ = 1п|| + Е|-; Е0 = е1; Е = е1 (|-1); Еп = £пе£- пКп_,, п = 1,2,...; I п=1 п!

ге1 соъа£ Сп те1 вта| А ^ пРи а>1 ]---ё£ = Сг (а|) + ^ п+1 п ; ]---= ^г(а£) + ^

| ^ «п+1(п + 1)Г | £0 ап+\п +1)!'

С0 = ьта!; Сп = ап£п Бта|-п£п-1;

=-соваЦ; £я = -ап£п оо$,а! + пСп_х, п = 1,2,....

В окончательном виде точное решение линеаризованной задачи имеет вид:

V и х = х0 соб«t + вт«0t +--(а\1\ + а212 - ЪХ3Х - );

«0 «0

Vx = -х0«0 ^ Ф + К0 С08 «0г - и Ъ1 + Ъ212 + а1З1 + а232 ) ;

у = у0ск«0 +-->?кю0лр21 -1)--= (а313 -а^^ -Ъъ33 -ЪЗ ^);

«0л/2 2«0 « 2У12

у =«0 у^2"к«0»21 1 'у0"к«0м2 - и (а313 - а41 4 - Ъ3З 3 Ъ4З 4

а

V = «оУо42 + Ку^ка0421--sка042t - и(а313 -- ЪъЗъ - ЪЗ4) .

0*2 2

Это решение характеризуется ошибками ~10 км по местоположению и ~120 м/с по скорости при снижении быстродействия в несколько раз по сравнению с методами [1,2] в зависимости от числа удерживаемых членов разложения подынтегральных выражений.

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

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

Рассмотрим теперь на этом же примере численные методы с учетом быстродействия их программной реализации. Предварительно отметим, что время выполнения компьютерных расчетов существенно зависит от алгоритмических и программных приемов повышения скорости вычислений [4]. Так, например, операция возведения числа во вторую степень выполняется на два порядка медленнее, чем умножение числа само на себя. Только одно это удваивает время выполнения шага интегрирования методом Рунге-Кутты, если при компиляции программного кода не используется режим оптимизации.

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

функций. При интегрировании уравнений движения с постоянным шагом h для случая линейной программы угла тангажа число вычислений этих функций сокращается до минимума, если предварительно рассчитать их значения в начальный момент времени, а затем каждое последующее определять через предыдущее по известным формулам тригонометрии:

sin (+1 = sin(( + (h /2) = sin ( cos((h /2) + cos ( sin((h/ 2);

cos (+1 = cos(( + (h/2) = cos ( cos((h /2) - sin ( sin((h / 2).

В целом, численное решение уравнений движения методом Рунге-Кутты, реализованное в виде специализированной, ориентированной на максимальное быстродействие подпрограммы, сокращает время выполнения примерно в 3 раза.

На внеатмосферном участке уравнения движения не зависят от скорости полета. Поэтому здесь более эффективны методы интегрирования дифференциальных уравнений второго порядка, не зависящих от первой производной. Например, трехэтапный метод Нюстрёма [5] при меньшем объеме и времени вычислений имеет такой же порядок точности как метод Рунге-Кутты. А при равном с методом Рунге-Кутты быстродействии в 1,5 раза более высокой точностью обладает двухэтапный метод, вычислительная схема которого может быть получена обобщением на случай дифференциальных уравнений второго порядка приведенного в [5] метода для уравнений первого порядка, с той особенностью, что первый расчет правых частей дифференциальных уравнений производится в точке, смещенной относительно начала шага. В соответствии с этим методом кинематические параметры траектории R = (x, y)T и V = (Vx ,V )т в конце шага интегрирования определяются

соотношениями:

V = V° + (q + q-)/6;

R = R0 + hV0 + h(a2 q + ax q2) / 4, где q = hf (a h, R + ax h V0);

q = hf (a2h, R + a2hV0 + hq /3);

a= (1-1/V3)/2; a2= (1 + 1/V3)/2;

f (t, R) - вектор-функция правых частей системы уравнений движения.

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

V = V + Ж + (£ + 4^2 + /6; Я = Я + £ + ЬУ, + (Ё! + 2 Ё2)й2/6;

Ё1 = Ё (Яо);

Ё2 = Ё[Яо + VоЬ /2 + Ё1Ь2 / 8 + £(Ь / 2)];

Ёз = Ё[Яо + VоЬ + Ё2Ь2 /2 + £(Ь)]. В этих соотношениях гравитационное ускорение Ё = (£х, ё )Т рассчитывается в начальной, промежуточной и конечной точках шага интегрирования, а вектора кажущейся скорости Ж = (Жх,Жу )Т и кажущегося пути £ = (£х, £ )Т определяются как интегралы от проекций кажущегося ускорения:

г г г

Ж = | Ж (г)йг = Г йг; £х = | Жх (г)йг;

о о о

г г р • г

Жу = { Ж (г )йг ={ -81пр йг; £у = \ажу (г )йг.

В рассматриваемой задаче они выражаются аналитически:

Р со&(<р0+фг) С рТ рТ ^

т - тг

ТЖГ Г-СОК(ро +р1) ^ .

Жх =1- .-йг = и соъртСг\ + Б1пр

1 т ^ р(т-г) р(т-г) у

о '"о

ш \Р81п(Ро +рг)^ ( . г.рр р ^

Ж =1-о-йг = и ътртСгп -соърт

т0 - тг

\ Р(Т-г) Р(Т-г)у

Б1п(р + рг) - б1п р0

(1)

= а — - Жх (Т - г);

Р

£ = и СОВ(ро +рг) - СОБ ро - Ж (Т - г).

у р у

Для вырожденного случая р = о эти формулы принимают вид:

Т

Жх = Ж соБр0; Жу = Ж Б1пр0; Ж = и 1п-;

Т — г

£х = £ соБр; £ = £ Б1пр0; £ = иг - Ж(Т - г).

По сравнению с методом Рунге-Кутты полуаналитический модифицированный трехэтапный метод Нюстрёма на порядок улучшает точность определения координат и на два порядка - вектора скорости, но занимает почти в 8 раз больше процессорного времени. Используя вместо обращения к стандартным подпрограммам вычисления интегральных синусов и косинусов следующие соотношения:

о

о

о

о

¿£ =Я£2 - =Уз(£22п-1 -£2п-1), з=—(-1п^1— I £ £ £ £ £ '(£ £ ^ ' (2п -1)! (2п -1)

^-а£ = 1п£(£22п-£2п), с=

£ £1 £1 (2п)!2п

можно сократить затраты времени на 40%. (Здесь коэффициенты при степенях переменной вычисляются одновременно для верхнего и нижнего пределов интегрирования, и требуется только одно «медленное» вычисление логарифма).

Еще более существенно (в 4 раза) повышается быстродействие, если, исходя из требуемой точности, заранее рассчитать необходимое число коэффициентов з и с, а суммы степенных рядов вычислять по схеме Горнера:

N

£з1£2п-1 = £(3 + £2(з2 + £2(з + £2^ ... + £2sN).);

п=1

N

£с£2п = £2(с, +£2(с2 + £2(съ +£2(сА... + £2CN).).

п=1

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

Е = -4 * Я3

Е = - еУ - 3УКЁ);

Я

Е = -Е(ЯУ - 6У,У ) + 3(5К,2 - V2 - Яа,)Е],

Я

где 4 - гравитационный параметр Земли;

2 = ,/7^)2; V, = ; „ =

Я Я

Представив решение уравнений движения как

г

V = V + Г + Ге , ГЕ =| Е (2 (т;

0

г

Я = Я0 + у0г + ^ + Яе , ЯЕ = \У^ ,

и используя первые производные ^, , , вычисленные в точках:

Я - Ц,; V - Г0;

- - -к - к2 - к3 - - — к ■ к2

Я = Я + V, - + £ (к/2) + £ — + £ —; V = V, + Ж (к/2) + £ - + £ —;

- - - - к2 • • к3 - - — ■■ к2

Я - Я, + Voк + £ (к) + (3£ + 2£2) — + (£ + £)—; V = V, + Ж (к) + £к + (£ + 2£) —,

10 30 6

можно найти уточненные гравитационные составляющие векторов скорости и положения:

к

V, - — [14Ё! + 32Ё2 +14§з + к(Ё! -£)];

- .1,Л1 , 32Л9 + 14е, +

6 601 1 2 3

к2 (3)

Я, - £ +112£2 +19§3 + к(5£ - В|2 - 2,3)].

^ ^0['9 61 1126 2 + 19 £ 3 + к(561 8§ 2 3.

Это позволяет уменьшить погрешность координат с 1052 м до 19 м и скорости с 0,888 м/с до 0,324 м/с.

Очевидно, что повышение точности расчетов сопровождается снижением быстродействия. В данном случае время вычислений увеличивается на 30%. Вдвое меньшие дополнительные затраты времени требуются для метода типа «прогноз-коррекция», в котором учитываются вторые производные от гравитационного ускорения.

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

К и Я :

-тт* ___:_ t .11 t — * _ t _:_ t .11 t 4

V = got + 60 — + 60-; Я = §0- + 60- + ,0— t.

После этого вычисленные по формулам (2) кинематические параметры траектории первого приближения позволяют рассчитать гравитационное ускорение § * и его производные § *, в найденной прогнозируемой конечной точке траектории. Построив степенной полином для § ^) с этими параметрами как граничными условиями, и проинтегрировав его, запишем окончательные выражения для скорректированных величин V и Я^ :

V - 1 6 2

- t2 Я - — 6 14

Ё0 + § + "1 §0 - § + — (6"0 + §)

5§0 + 2§* I 13§0 -8§* + (4+ 3§*)

(4)

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

логарифмическом масштабе на рисунках 1 и 2. На каждой кривой, относящейся к какому либо методу, крайняя левая точка соответствует решению уравнений движения за один шаг интегрирования, а последующие точки - численному интегрированию за несколько шагов. Программная реализация каждого метода выполнена с применением приемов повышения скорости вычислений. Время расчета дано в процентах по отношению к методу Рунге-Кутты четвертого порядка точности с шагом интегрирования И =1 с, а полученные с его помощью результаты приняты за эталонные.

^ 4 § 104

О О И ей

н и

^ „ „2 О 10

то

Л Н О О К

В

<и 0

Ер10

о С

• о о приближенный метод из работы [1] ° приближенный метод из работы [2] -О— трехэтапный метод Нюстрёма -О— смещенная начальная точка —□— модифицированный метод Нюстрёма —□— метод "прогноз-коррекция" —□— аппроксимация гравитационного ускорения —Д— аппроксимация + расчет производных

0.5

1.5

время расчета, %

2.5

Рисунок 1. Зависимость погрешности расчета координат от затрат процессорного времени

0

1

2

3

о

и

К 10:

н

о

о

О

И

о

ей

н

и сг 10(

о

ей

Л

Н

О

о

К

3

и Ер 10"

о

с

о о приближенный метод из работы [1] ° приближенный метод из работы [2] -О— трехэтапный метод Нюстрёма —О— смещенная начальная точка —□— модифицированный метод Нюстрёма —□— метод "прогноз-коррекция" —□— аппроксимация гравитационного ускорения —Д— аппроксимация + расчет производных

0.5

1.5

время расчета, %

2.5

Рисунок 2. Зависимость погрешности расчета скорости от затрат процессорного времени

2

0

0

1

2

3

Из приведенных данных видно, что приближенное (за 1 шаг) численное интегрирование уравнений движения трехэтапным методом Нюстрёма и методом со смещенной начальной точкой обеспечивает более высокую точность решения по сравнению с приближенными линеаризованными методами, не уступая им в быстродействии. Модифицированный метод Нюстрёма и метод «прогноз-коррекция» повышают точность расчетов на два порядка, увеличивая затраты процессорного времени в 1,5-2 раза. По критерию «точность-быстродействие» лучшие результаты при повышенных требованиях к точности имеет последний из них. В представляющем практический интерес диапазоне погрешностей расчета скорости 0,01-1 м/с и 1-100 м по координатам вычислительная сложность рассмотренной задачи примерно одинакова у всех рассмотренных численных методов, то есть при равных требованиях к точности затраты процессорного времени у них отличаются не более чем в 2 раза. Причем в отношении быстродействия не меньшее значение чем математическая структура конкретного метода имеет качество ее программной реализации, то есть оптимизация схемы вычислений и использование алгоритмических и программных способов повышения быстродействия. В целом, численное решение уравнений движения на внеатмосферном активном участке траектории допускает большой шаг интегрирования порядка 100 и более секунд, обеспечивая приемлемую для задач проектной баллистики точность и не ограничивая диапазон возможных траекторий.

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

Для полуаналитического модифицированного метода Нюстрёма производные от векторов скорости и положения определяются следующими формулами:

7' = Щ + (4^ + Ц= Щ+ 2^2/6; й = ф0,ф. (5)

Здесь производные от кажущихся параметров получаются дифференцированием формул (1):

Ж' = (-Ж )т; Г= (-£ , Бх )т;

п

Ж' = ' + Ж г

- Ж'(т -г) + -

и

(р

иг2

- Ж(т - г) + —

фг ео8(п + фг) - §т(п0 +фг) + Б1П п фг б1П(п + фг) + СОБ(П + фг) - СОБ п

при ф Ф 0;

СОБП

при ф = 0.

<?Х,Х у Я2 - 3х2 - 3ху

х , у Я5 - 3ху Я2 - 3 у2

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

Е'гл= Ё'г&щ (Н/2); ^ = ^ (ё2^2/2 + Щ (Н)), ' = Ро,Ф;

Як =

Вычисление производных в полтора раза увеличивает общие затраты процессорного времени, но это в два раза быстрее трехкратного расчета траектории этим же методом при определении производных разностным способом. При этом погрешность расчета производных не превышает 2%. Вдвое точнее находятся производные методом «прогноз-коррекция», но вследствие значительно более высокой сложности выражений для производных от формул (4) время расчета более чем удваивается по сравнению с модифицированным методом Нюстрёма, у которого за два шага интегрирования погрешность производных снижается до 0,2%.

Учитывая, что матрицу ^ можно использовать и для определения производных по

времени £ = ^Я, целесообразно одновременно с расчетом производных по параметрам программы угла тангажа п и ф применять формулы (3) для уточнения гравитационных составляющих Я и V , причем практически без дополнительного снижения

<

быстродействия и без необходимости уточнения производных от самих Я и V , так как

достигнутый в модифицированном методе Нюстрёма уровень погрешности формул (5) вполне достаточен для решения краевой баллистической задачи.

Соответствующие показатели «точность-быстродействие» этого метода вместе с характеристиками ранее рассмотренных показаны на рис.1 и рис.2. Мало уступая по скорости вычислений, но обеспечивая при этом расчет производных от кинематических параметров траектории по параметрам управления, этот метод сокращает время решения краевой баллистической задачи не менее чем в 2-3 раза по сравнению с другими рассмотренными методами, и показывает наиболее высокую точность решения за один шаг интегрирования.

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

1. Аппазов Р.Ф., Сытин О.Г. Методы проектирования траекторий носителей и спутников Земли. - М.: «Наука», 1987. - 440 с.

2. Баллистические ракеты и ракеты-носители. / Под ред. О.М.Алифанова. - М.: «Дрофа», 2004. -512 с.

3. Инженерный справочник по космической технике. /Под общей редакцией А.В.Соловьева. М.: «Воениздат», 1969. -696 с.

4. Немнюгин М.А., Стесик О.Л. Современный Фортран. СПб.: БХВ-Петербург, 2005. -496 с.: ил.

5. Хайрер Э., Нёрсетт С, Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи: Пер. с англ. — М.: Мир, 1990. —512 с, ил.

6. Беневольский С.В., Козлов П.Г. Полуаналитический метод восстановления траекторий ЛА по обобщенным проектным параметрам и параметрам системы управления и перспективы его использования. //Электронное научно-техническое издание «Наука и образование», №10, октябрь 2011г., http: \technomag .edu.ru/doc/216895. html

7. Беневольский С.В., Гончаренко В.И. Перспективные технологии решения краевых задач доставки отделяющихся частей ракет-носителей с терминальным наведением. //Авиационная и ракетно-космическая техника. Вестник Самарского государственного аэрокосмического университета. № 2 (26) 2011 г., с. 18-28.