Какова функциональная форма r(t) для спуска солнечного паруса с орбиты на Солнце?

Этот ответ предполагает, что космический корабль все время остается на почти круговой орбите. У нас есть

$$\frac{dE}{dt} = F_\tau v,$$ куда $E$энергия космического корабля (потенциальная + кинетическая), $v$это его скорость и $F_\tau$– тангенциальная составляющая силы солнечного давления. Для круговой орбиты $E\propto -\frac1r$(так $\frac{dE}{dt}\propto \frac{1}{r^2}\frac{dr}{dt}$) и $v\propto\frac{1}{\sqrt{r}}$. Пока угол между парусом и направлением на солнце остается постоянным (независимо от того, оптимален он или нет), $F_\tau\propto -\frac{1}{r^2}$. Итак, у нас есть $$\frac{1}{r^2}\frac{dr}{dt} \propto -\frac{1}{r^{5/2}},$$ или $$\frac{dr}{dt} \propto - \frac{1}{\sqrt{r}}.$$

Решение такого дифференциального уравнения имеет вид

$$r(t) \propto (T - t)^{2/3},$$ куда $T$(что больше времени пуска) определяется условиями пуска и коэффициентом пропорциональности в уравнении.

Похоже, ответ @Litho удался !

$$r(t) \propto (T - t)^{2/3}$$

Я провел быструю симуляцию, основанную на спецификации LightSail 2 Планетарного общества: 5 кг и площадь паруса 32 м^2. Я установил его под углом отражателя 45 °, так что давление солнечного света приводит к радиальной силе наружу (импульс падающего света) плюс тангенциальная сила, направленная вперед (импульс отраженного света).

Ускорение, вызванное импульсом света, падающего на поверхность или покидающего ее, равно

$$\frac{AI}{mc} = \frac{AI_0}{mc} \left(\frac{\text{1 AU}}{r}\right)^2$$

куда$I_0$- солнечная постоянная (интенсивность на 1 а.е.) около 1361 Вт/м^2. Смотрите этот ответ , чтобы узнать больше о солнечном давлении и ускорении солнечными парусами. Не забудьте разделить площадь паруса на$\sqrt{2}$чтобы получить площадь проекции под углом 45°.

Я начал с круговой орбиты на расстоянии 1 а.е. и проинтегрировал 15,35 года.

Оказывается, что$T$время прибытия, поэтому на первом графике я просто сравниваю$r$, расстояние до Солнца в моделировании, к простому выражению:

$$\text{1 AU} \left(1 - \frac{t}{T}\right)^{2/3}$$

и вуаля идеальная посадка! Покачивания происходят из-за того, что я начал с гелиоцентрической круговой орбиты в 1 а.е. и скорости$\sqrt{GM_{Sun}/1 AU} =$29783 м/с с эффектами солнечного давления в полную силу (замедление, небольшая внешняя сила, уменьшающая гравитацию, поэтому орбита очень слегка эллиптическая.

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

$$-\frac{GM}{r^2} \ \ \text{and} \ \ +\frac{AI_0}{\sqrt{2}mc} \frac{\text{1 AU}^2}{r^2}$$

Численно на расстоянии 1 а.$1/r^2$отношение двух фиксировано и не зависит от расстояния. В этом случае соотношение составляет около 289:1.

def deriv (X, t):

    r,  v  = X.reshape(2, -1)
    nr, nv = [thing / np.sqrt((thing**2).sum()) for thing in (r, v)] # normals
    rsqAU  = (r**2).sum() / AUsq

    acc_g     = -GMs * r * ((r**2).sum())**-1.5
    acc_solar = (Area/np.sqrt(2.) * I_zero / (m * c) / rsqAU) * (nr - nv) # radially out, and prograde

    return np.hstack((v, acc_g + acc_solar))

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

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
degs, rads        = 180/pi, pi/180

AU     = 1.495978707E+11       # m
AUsq   = AU**2                 # m^2

GMs  = 1.327E+20               # m^3/s^2

km   = 1000.                   # meters
year = 365.2564 * 24. * 3600.  # seconds

# http://www.planetary.org/explore/projects/lightsail-solar-sailing/lightsail-faqs.html
m      = 5.                    # kg
c      = 3E+08                 # m/s
I_zero = 1361.                 # 1361 W/m^2 (at 1 AU)
Area   = 32.                   # m^2

time = np.arange(0, 15.35*year, 1E+05)  # seconds

v0    = np.sqrt(GMs/AU)

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

print "X0: ", X0

answer, info = ODEint(deriv, X0, time, rtol=1E-10, full_output=True)

print answer.shape

x, v = answer.T.reshape(2, 2, -1)
r    = np.sqrt((x**2).sum(axis=0))
x, y = x

if True:
    plt.figure()
    plt.subplot(2, 1, 1)
    plt.plot(x/km, y/km)
    plt.title('heliocentric de-orbit (km)')
    plt.subplot(2, 1, 2)
    plt.plot(time/year, x/km)
    plt.plot(time/year, y/km)
    plt.plot(time/year, r/km, '-r', linewidth=2)
    plt.title('x, y and r (km) vs time (years)')
    plt.show()

if True:
    T0 = time.max()
    plt.figure()
    plt.plot(time/year, r/km)
    plt.plot(time/year, AU*(1-time/T0)**(2./3)/km)
    plt.title('r and  AU(1-t/15.35)^(2/3) (km) vs time (years)')
    plt.show()