Уравнение для скорости и расстояния от Солнца космического корабля с солнечным парусом

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

Даже в самых простых случаях вы говорите о разложении эллиптических интегралов в ряды или непрерывные дроби. Однако, если у вас есть книга Баттина по астродинамике, он показывает, как делать это вручную, как это делали такие люди, как Гаусс и Лежандр, и что должны были пройти все его студенты в Массачусетском технологическом институте, на страницах 408–418. Тем не менее, он начинает раздел с

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

Обзор еще нескольких методов, учитывающих, как именно вы планируете двигаться по спирали, которые я считаю в целом более доступными и полезными в данном случае, — это Петропулос и Симс, 2002 г. , «Обзор некоторых точных решений плоских уравнений». движения толкающего космического корабля». Обратите внимание, что они больше всего заинтересованы в использовании солнечных батарей для питания ионного двигателя, позволяющего вам двигаться в наиболее выгодном направлении, которое очень далеко от радиального!

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

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

Математика

Учитывая начальное ускорение на 1 а.е.$8.17 \times 10^{-4} m/s^2$.

Наклоните его на 45 градусов, чтобы сделать тягу тангенциальной, разделите на$\sqrt{2}$поскольку теперь он наклонен к Солнцу и учитывает падение с расстоянием от Солнца:

$$a_0 = 5.78 \times 10^{-4} m/s^2$$

$$a(r) = a_0 \frac{AU^2}{r^2}$$

в прямом направлении (в том же направлении, что и скорость течения).

уравнение vis-viva :

$$v^2 = \frac{GM}{r}$$

Сейчас

$$\frac{dv}{dt} = -\frac{a_0 AU^2}{GM^2} v^4$$

где$GM$стандартный гравитационный параметр Солнца $1.327 \times 10^{20} m^3/s^2$. Знак минус появляется, потому что мы знаем, что вопреки первому инстинкту, когда у нас есть ускоряющая сила в прямом направлении, мы, вопреки интуиции, замедляемся на ту же величину. Это также цитируется в нескольких других сообщениях здесь, я поищу другие ответы, чтобы процитировать...

Перепишите и решите:

$$\frac{dt}{dv} = -\frac{GM^2}{a_0 AU^2} v^{-4}$$

$$t(v) = t_0 - \frac{GM^2}{4 a_0 AU^2} v^{-3}$$

Если мы установим$t(v_0) = 0$другими словами, в нулевой момент времени мы движемся с орбитальной скоростью$v_0 = \sqrt{GM/AU}$ получаем 0,408 года.

Численно

Интересно, что когда я пытаюсь смоделировать то же самое численно, я получаю время выхода 0,515 года! Меня это удивляет по двум причинам:

1. Скорость никогда не достигает нуля, а это означает, что приведенная выше математика неверна!
2. И все же это довольно близко!

Заключение

На побег уйдет примерно полгода, но полного уравнения у меня пока нет.

Дальнейшее чтение:

• Можно ли достичь Солнца, не расходуя топливо/реакционную массу?
• Каков оптимальный угол схода солнечного паруса с орбиты к Солнцу с учетом радиальной тяги? (в настоящее время без ответа)

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

AU = 150E+06 * 1000   # meters
GM = 1.327E+20   # m^3/s^2
a0 = 8.17E-04 / np.sqrt(2)  # dv/dt at 1 AU
year = 365.2564 * 24 * 3600 

coef = ( GM**2 / (4 * a0 * AU**2) )

v0 = np.sqrt(GM/AU)
print('v0: ', v0)
print('coef: ', coef)

t0 = coef / v0**3
print('t0 / year: ', t0 / year)

def deriv(X, t):
    x, v = X.reshape(2, -1)
    vhat = v / np.sqrt((v**2).sum())
    rsq = (x**2).sum()
    acc_thrust = vhat * a0 * rsq / AU**2
    acc = acc_thrust - GM * x * rsq**-1.5
    return np.hstack((v, acc))

X0 = np.array([AU, 0, 0, v0])

time = np.linspace(0, 0.6, 10001) * year  # half year

answer, info = ODEint(deriv, X0, time, full_output=True)

print(answer.shape)

x, y, vx, vy = answer.T
r, speed = [np.sqrt((thing**2).sum(axis=0))
            for thing in answer.T.reshape(2, 2, -1)]

E = 0.5 * speed**2 - GM/r
E_norm = np.abs(E[0])

i_esc = np.argmax(E>=0)
things = time, r, x, y, speed, E
t_esc, r_esc, x_esc, y_esc, s_esc, E_esc = [thing[i_esc]
                                            for thing in things]

print('t_esc / year: ', t_esc / year)

if True:
    plt.figure()
    plt.subplot(2, 2, 1)
    plt.plot(time/year, r/AU)
    plt.plot([t_esc/year], [r_esc/AU], 'ok')
    plt.ylabel('r/AU')
    plt.xlabel('time (years)')
    
    plt.subplot(2, 2, 2)
    plt.plot(time/year, speed/1000)
    plt.plot([t_esc/year], [s_esc/1000], 'ok')
    plt.ylabel('speed (km/s)')
    plt.xlabel('time (years)')

    plt.subplot(2, 2, 3)
    plt.plot(time/year, E/E_norm)
    plt.plot([t_esc/year], [E_esc/E_norm], 'ok')
    plt.plot(time/year, np.zeros_like(time), '-k')
    plt.ylabel('Energy (norm)')
    plt.xlabel('time (years)')

    plt.subplot(2, 2, 4)
    plt.plot(x/AU, y/AU)
    plt.plot([x_esc/AU], [y_esc/AU], 'ok')
    plt.plot([0], [0], 'oy')
    th = np.linspace(0, 2*np.pi, 201)
    plt.plot(np.cos(th), np.sin(th), '-r', linewidth=0.5)
    plt.ylim(-1, 1.5)
    plt.gca().set_aspect('equal')
    plt.xlabel('AU')
    plt.ylabel('AU')
    plt.show()