4
$\begingroup$

This answer describes a spacecraft in a heliocentric orbit using a solar sail to "deorbit" into the Sun by reflecting sunlight into roughly the prograde direction. In comments it was proposed that while the solar flux and therefore thrust of the sail increases as the craft gets closer to the Sun, the delta-v required for a given change in orbit size increases as well.

For a spacecraft of mass $m$ and perfectly reflecting sail area $A$, what does the plot of orbit radius versus time $r(t)$ look like? Is it a straight line, power law ($p \neq 1$), exponential, logarithmic, or something else?

Since the answer mentions that a reflector angle of 45°, resulting in exactly prograde thrust is not necessarily the optimum, you are welcome to consider either 45° or the optimum angle, fixed or dynamic (though I don't think it would change with distance).

"bonus points" for the derivation of $r(t)$!

$\endgroup$
  • 2
    $\begingroup$ Well-- if a craft like this is ever designed it absolutely must be called ICARUS. $\endgroup$ – Magic Octopus Urn Aug 25 '18 at 4:16
  • $\begingroup$ Probably solution is not an analytical function. $\endgroup$ – Heopps Aug 25 '18 at 6:10
  • 1
    $\begingroup$ @Heopps for a circular orbit or a slow spiral, I'll bet it is! Let's see what happens... $\endgroup$ – uhoh Aug 25 '18 at 7:17
2
+100
$\begingroup$

This answer assumes that the spacecraft stays in an almost-circular orbit all the time. We have $$ \frac{dE}{dt} = F_\tau v, $$ where $E$ is spacecraft's energy (potential + kinetic), $v$ is its speed, and $F_\tau$ is the tangential component of the solar pressure force. For a circular orbit, $E\propto -\frac1r$ (so $\frac{dE}{dt}\propto \frac{1}{r^2}\frac{dr}{dt}$) and $v\propto\frac{1}{\sqrt{r}}$. As long as the angle between the sail and the direction to the sun stays constant (whether it's optimal or not), $F_\tau\propto -\frac{1}{r^2}$. So we have $$ \frac{1}{r^2}\frac{dr}{dt} \propto -\frac{1}{r^{5/2}}, $$ or $$ \frac{dr}{dt} \propto - \frac{1}{\sqrt{r}}. $$

A solution of a differential equation like this has the form $$ r(t) \propto (T - t)^{2/3}, $$ where $T$ (which is greater than the starting time) is determined by the starting conditions and the proportionality coefficient in the equation.

$\endgroup$
  • $\begingroup$ okay it looks like the onus is on me to check it, I'll give it a go in the morning. Thanks! $\endgroup$ – uhoh Aug 27 '18 at 19:16
  • $\begingroup$ @uhoh Sorry, I made a mistake. The power should be 2/3, not 2. $\endgroup$ – Litho Aug 28 '18 at 2:49
  • $\begingroup$ That looks better! It also appears to have the same shape as my simulation (I'll post something soon). Also, I think $T$ is "greater than (or equal to) the finish time" as well. Exponentiation to fractional powers returns complex numbers for negative arguments, for example $-1^{2/3}=-1/2+\sqrt{3/4}j$ So probably $T$ may end up being just the arrival time itself, since $r(t)=0$. $\endgroup$ – uhoh Aug 28 '18 at 4:00
0
$\begingroup$

Looks like @Litho's answer nailed it!

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

I did a quick simulation based on the Planetary Society's LightSail 2's spec of 5kg and 32 m^2 sail area. I fixed it at a reflector angle of 45° so that the pressure from sunlight results in a radial force outward (momentum of the incident light) plus a a tangential force prograde (momentum of the reflected light).

The acceleration due to the momentum of light either striking or leaving a surface is just

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

where $I_0$ is the solar constant (intensity at 1 AU) of about 1361 W/m^2. See this answer for more on solar pressure and acceleration by solar sails. Remember to divide the area of the sail by $\sqrt{2}$ to get the projected area at 45°.

I started in a circular orbit at 1 AU and integrated for 15.35 years.

It turns out that $T$ is the arrival time, so in the first plot I just compare $r$, the distance to the Sun in the simulation, to the simple expression:

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

and voilà a perfect fit! The wiggles are due to the fact that I started with a heliocentric circular orbit of 1 AU and a velocity of $\sqrt{GM_{Sun}/1 AU} =$ 29783 m/s with the solar pressure effects at full strength (deceleration, slight outward force reducing gravity, and so the orbit is very slightly elliptical.

Radial accelerations due to the Sun's gravity and incident radiation pressure are given by:

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

Numerically at 1 AU they are 5.930E-03 and 2.053E-05 m/s^2 respectively, and because both scale as $1/r^2$ the ratio of the two is fixed and independent of distance. In this case the ratio is about 289:1.

enter image description here

enter image description here

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()
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.