2. Программа lrkla 1

Для расчета параметров орбиты перелета ρ, e и составлена программа LRKLA 1.

Перед запуском программы необходимо ввести исходные данные:

RФ= - длина начального радиуса-вектора, км;

RK= – длина конечного радиуса-вектора, км;

FF=Ф – угловая дальность, рад;

Q=μ – гравитационный параметр, /;

TPOL= - время перелета, сутки;

EPS= β – погрешность вычисления полуоси , км.

В программе реализован рассмотренный выше алгоритм.

Программа позволяет производить расчёт для трёх типов орбит: эллиптической, гиперболической и параболической.

В результате выполнения программы вычисляются параметры орбиты: p - фокальный параметр; E=e - эксцентриситет; A=a -большую полуось(км); а также – скорость в начальной точке; -скорость в конечной точке.

Скорость для начальной и конечной точки орбиты вычисляется из закона сохранения энергии:

Блок-схема программы LRKLA1 приведена на рис. 1.3.

В программе для решения трансцендентных уравнений(1.1), (1.2) методом хорд используется подпрограмма TRANS, из которой происходит обращение к подпрограммам GRAN и YRAN.

Подпрограмма GRAN служит для задания минимального a₁ и максимального значения a₂ полуоси орбиты по формулам (1.6), (1.7).

Подпрограмма YRAN вычисляет функцию невязки y(a) (см. формулу(1.8)).

Параметры орбиты вычисляются по формулам (1.3)- (1.5) с помощью подпрограммы PARAM, а также скорости и по формуле (1.10).

3. Порядок выполнения работы.

Требуется рассчитать траекторию перелёта с круговой орбиты Земли на круговую орбиту Марса. Будем считать, что орбиты планет и траектория перелёта лежат в одной плоскости. Орбиты Земли и Марса считаем круговыми.

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

На этапе движения КА в грависфере Земли решается задача двух тел Земля- КА. Этап заканчивается в момент выхода КА из грависферы Земли. Этот участок траектории называется геоцентрическим.

После выхода из грависферы Земли движение КА происходит в грависфере Солнца. На этом участке, называемом гелиоцентрическим, траектория КА исследуется в рамках задачи двух тел Солнце- КА. Начальными условиями для гелиоцентрического участка являются характеристики траектории в конце геоцентрического участка.

После выхода траектории в грависферу Марса решается задача двух тел Марс- КА.

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

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

При расчёте гелиоцентрического участка полёта удаётся решить основную задачу полёта- попасть в окрестность планеты назначения.

Расчёт гелиоцентрического участка полёта производится в предположении, что радиусы гравитационных сфер планет равны нулю. В этом случае считается, что межпланетный участок траектории КА начинается в точке, где располагается центр Земли (а самой Земли и её грависферы нет). Гелиоцентрический участок траектории заканчивается в точке, где располагается планета назначения- Марс.

Время начала гелиоцентрического участка траектории считается равным времени старта КА с орбиты искусственного спутника Земли. Аналогичные допущения делаются при анализе времени полёта КА к планете цели.

Рассмотрим расчёт гелиоцентрического участка траектории (рис 1.4).

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

Время перелета с орбиты 1 на орбиту 2 в этом случае определяется формулой

= π (1.11)

Для встречи с Марсом в точке необходимо. Чтобы в момент старта Марс находился в точке , смещенным относительно точки встречи на угол = (2 π/)∙57,3, где = 687 сут. – период обращения Марса вокруг Солнца.

Учитывая, что =132,718∙ /; = =149,6∙ км;

= = 227,9∙ км, из (1.11) получаем = 259 сут.

За время полета КА по гомановской траектории Земля сместится на угол

= (2 π/)∙57,3 = , где = 365 сут.

Так как реализация оптимальных полетов требует вполне определенного расположения Марса относительно Земли, а такие ситуации возникают довольно редко (через 2,14 года), то перелеты происходят не по энергетически оптимальным траекториям.

При выполнении лабораторной работы нужно определить параметры траектории по данным, приведенным в таблице 1.1, а также параметры ρ и e оптимальной гомановской траектории. Кроме того, необходимо исследовать влияние времени перелета на энергетику старта с орбиты спутника Земли и перехода на орбиту спутника Марса.

Таблица 1.1

Номер варианта	Ф, рад	, сут.	, сут.	, сут.
4	2,3	180	230	280

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

tп, сут	V0, км/с	Vk, км/с	p	e	a, км
180	33.0340	21.9463	1.7413e+008	0.3222	1.9430e+008
230	32.2767	20.7891	1.3526e+008	0.5034	1.8118e+008
280	30.1416	19.4361	1.2047е+008	0.7080	1.3301e+008

Рисунок 1. Траектория при Тп = 180 сут

Рисунок 2 Траектория при Тп = 230 сут

Рисунок 3. Траектория при Тп = 280 сут

Код программы:

function mars()

% ввод начальных значений

r0 = 149.6*10^6; %длина начального РВ, км

rk = 227.9*10^6; %длина конечного РВ, км

Fi = 2.3; %угловая дальность, рад

mu = 132.718*10^9; %грав. пар-р Солнца

tp1 = 280*3600*24; %время перелета, с

eps = 8; %погрешность вычисления большой полоуоси, км

i=0;

S = sqrt(r0^2+rk^2-2*r0*rk*cos(Fi));

tpar = ((r0+rk+S)^(1.5)-((r0+rk-S)^(1.5))*sign(sin(Fi)))/(6*sqrt(mu));

del = (tp1-tpar)/3600/24;

if(del<0) %орбита гиперболическая

a1 = r0/50; %минимальное значение большой полоуоси

a2 = rk*10; %максимальное значение большой полоуос

ya3 = 9000;

while(ya3>6*10^(-9))

a3 = a1 - (y(r0,rk,S,Fi,mu,a1,tp1)/(y(r0,rk,S,Fi,mu,a2,tp1)-y(r0,rk,S,Fi,mu,a1,tp1)))*(a2-a1);

ya3 = y(r0,rk,S,Fi,mu,a3,tp1)

if(ya3*y(r0,rk,S,Fi,mu,a1,tp1)<0)

a2=a3;

else

a1=a3;

end

end

%break;

else% орбита эллипт.

a1 = (r0+rk+S)/4;

a2 = 6*rk;

ya3 = 9000;

while(abs(a2-a1)>eps)

a3 = a1 - (y2(a1)/(y2(a2)-y2(a1)))*(a2-a1);

i = i+1;

ya3 = y2(a3);

if(y2(a3)*y2(a1)<0)

a2=a3;

else

a1=a3;

end

if(i>1000)

break;

end

end

a = a3

end

alpha = 2*asinh((r0+rk+S)/4*abs(a));

beta = 2*asinh((r0+rk-S)/4*abs(a));

gm = 2*asin(sqrt((r0+rk-S)/(r0+rk+S)));

tm = ((r0+rk+S)^(1.5))/(8*sqrt(mu))*(pi-sign(sin(Fi))*(gm-sin(gm)));

p = a/S^2*4*r0*rk*(sin(Fi/2))^2*(sin((alpha+sign(sin(Fi))*sign(tm-tp1)*beta)/2))^2

e = sqrt(1-p/a)

v0 = sqrt(mu*(2/r0-1/a))

vk = sqrt(mu*(2/rk-1/a))

%end

t=[0:pi/180:2*pi];

axis equal

plot(r0*cos(t),r0*sin(t),'k',rk*cos(t),rk*sin(t),'k');

hold on

plot(r0,0,'ro');

t = [0:pi/180:2*pi/687/24/3600*tp1];

plot(r0*cos(t)+(rk-r0)/1.64.*t.*cos(t),r0*sin(t)+(rk-r0)/1.64*t.*sin(t));

plot(rk*cos(2*pi/687/24/3600*tp1),rk*sin(2*pi/687/24/3600*tp1),'ro')

end

function F = y(r0,rk,S,Fi,mu,a,t)

alpha = 2*asinh((r0+rk+S)/4*abs(a));

beta = 2*asinh((r0+rk-S)/4*abs(a));

tp = sqrt((abs(alpha))^3/mu)*(sinh(alpha)-alpha-sign(sin(Fi))*(sinh(beta)-beta));

F = t-tp;

end

function F2 = y2(a)

r0 = 149.6*10^6; %длина начального РВ, км

rk = 227.9*10^6; %длина конечного РВ, км

Fi = 2.3; %угловая дальность, рад

mu = 132.718*10^9; %грав. пар-р Солнца

tp1 = 280*3600*24; %время перелета, с

S = sqrt(r0^2+rk^2-2*r0*rk*cos(Fi));

tpar = ((r0+rk+S)^(1.5)-((r0+rk-S)^(1.5))*sign(sin(Fi)))/(6*sqrt(mu));

gm = 2*asin(sqrt((r0+rk-S)/(r0+rk+S)));

g = 2*asin(sqrt((r0+rk-S)/(4*abs(a))));

e = 2*asin(sqrt((r0+rk+S)/(4*abs(a))));

tm = ((r0+rk+S)^(1.5))/(8*sqrt(mu))*(pi-sign(sin(Fi))*(gm-sin(gm)));

tp = ((a^(1.5))/sqrt(mu))*(pi+sign(tm-tp1)*(e-sin(e)-pi)-sign(sin(Fi))*(g-sin(g)));

F2 = tp1-tp;

end