¿Cuál es la forma funcional de r(t) para una vela solar que sale de órbita hacia el Sol?

Esta respuesta asume que la nave espacial permanece en una órbita casi circular todo el tiempo. Tenemos

$$\frac{dE}{dt} = F_\tau v,$$ donde $E$es la energía de la nave espacial (potencial + cinética), $v$es su velocidad y $F_\tau$es la componente tangencial de la fuerza de presión solar. Para una órbita circular, $E\propto -\frac1r$(asi que $\frac{dE}{dt}\propto \frac{1}{r^2}\frac{dr}{dt}$) y $v\propto\frac{1}{\sqrt{r}}$. Siempre que el ángulo entre la vela y la dirección al sol se mantenga constante (ya sea óptimo o no), $F_\tau\propto -\frac{1}{r^2}$. Entonces tenemos $$\frac{1}{r^2}\frac{dr}{dt} \propto -\frac{1}{r^{5/2}},$$ o $$\frac{dr}{dt} \propto - \frac{1}{\sqrt{r}}.$$

Una solución de una ecuación diferencial como esta tiene la forma

$$r(t) \propto (T - t)^{2/3},$$ donde $T$(que es mayor que el tiempo de inicio) está determinado por las condiciones de inicio y el coeficiente de proporcionalidad en la ecuación.

¡ Parece que la respuesta de @Litho lo logró !

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

Hice una simulación rápida basada en la especificación LightSail 2 de Planetary Society de 5 kg y 32 m ^ 2 de área de vela. Lo fijé en un ángulo de reflector de 45 ° para que la presión de la luz solar resulte en una fuerza radial hacia afuera (momento de la luz incidente) más una fuerza tangencial progresiva (momento de la luz reflejada).

La aceleración debida a la cantidad de movimiento de la luz, ya sea al incidir o al salir de una superficie, es

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

donde$I_0$es la constante solar (intensidad a 1 UA) de aproximadamente 1361 W/m^2. Consulte esta respuesta para obtener más información sobre la presión solar y la aceleración de las velas solares. Recuerda dividir el área de la vela por$\sqrt{2}$para obtener el área proyectada a 45°.

Empecé en una órbita circular a 1 UA y me integré durante 15,35 años.

Resulta que$T$es el tiempo de llegada, por lo que en la primera gráfica solo comparo$r$, la distancia al Sol en la simulación, a la expresión simple:

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

y voilà un ajuste perfecto! Los movimientos se deben al hecho de que comencé con una órbita circular heliocéntrica de 1 AU y una velocidad de$\sqrt{GM_{Sun}/1 AU} =$29783 m/s con los efectos de la presión solar a plena potencia (desaceleración, una ligera fuerza hacia el exterior que reduce la gravedad, por lo que la órbita es muy ligeramente elíptica.

Las aceleraciones radiales debidas a la gravedad del Sol y la presión de radiación incidente vienen dadas por:

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

Numéricamente a 1 UA son 5.930E-03 y 2.053E-05 m/s^2 respectivamente, y porque ambos escalan como$1/r^2$la relación de los dos es fija e independiente de la distancia. En este caso, la relación es de aproximadamente 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()