Equation for Velocity and distance from Sun of a solar sail powered spacecraft

Question

I have an interest in solar sails as a method of propulsion. I have thought of the following interesting question about the workings of solar sails, which I just don’t have the mathematical ability to be able to answer, but hopefully some of the members of this forum with a better mathematical knowledge may be able to assist me….

From https://en.wikipedia.org/wiki/Solar_sail, at the mean Earth-Sun distance (1 AU) radiation pressure from the Sun on a 90% efficient solar sails generates a thrust of $8.17$ $μN/m^2$ if the sail is perpendicular to the Sun’s direction. Suppose we have spacecraft in orbit at a distance of 1 AU from the Sun and let’s say the spacecraft weighs 10 kg. For simplicity, it is orbiting the Sun in a circular orbit. So, its orbital speed must be constant at ~29.8 km/s. The spacecraft then unfurls a solar sail, which is let’s say 10μ mm thick, 1000 square metres in area and made out of material of density $3000$ $kg/m^3$. So, the mass of the sail is only 0.03 kg, which is negligible compared to the mass of the spacecraft.

Therefore, the total outward force due to solar radiation pressure on the sail (assuming that it is normal to the direction of the Sun and it is 90% efficient) is $8.17$ x $10^{-3} N$. This gives an outward acceleration (force/mass) on the spacecraft of $8.17$ x $10^{-4} m/s^2$ away from the Sun. At a distance of 1 AU from the Sun the centripetal acceleration due to this Sun’s gravity $5.9$ x $10^{-3} m/s^2$. However, even though the net force on the spacecraft is still inwards towards the Sun, my intuition tells me the spacecraft will gradually spiral outward in its orbit around the Sun.

I would also assume that the strength of the outward force on spacecraft due to radiation pressure on the sails will (like the Sun’s gravity) fall as the inverse square of its distance from the Sun.

But the questions I am unable to answer are how do I calculate as a function of an arbitrary time t since the sail was unfurled.

• the spacecraft’s distance from the Sun ?
• And also its orbital velocity? (which I assume will be always normal to the Sun)

Answer 1

The math gets very complicated very quickly. Before even considering thrust, recall that the Kepler two-body orbit solution itself is not an explicit function of time: radius is given as a function of angle (true anomaly), but time is proportional to mean anomaly, which must be converted to eccentric anomaly and then true anomaly by numerical solution of Kepler's transcendental equation.

Even for the simplest cases, you're talking about series or continued fraction expansions of elliptic integrals. However, if you have Battin's book on astrodynamics, he shows how to do it by hand, the way people like Gauss and Legendre did, and which all his students at MIT were expected to fight through, on pages 408 to 418. That said, he begins the section with

One of the possible disturbing accelerations affecting a spacecraft is the thrust acceleration produced by the vehicle's engines. Trajectory determination under these circumstances, as in almost all problems of disturbed motion, generally requires the application of special numerical techniques. However, there are several examples of some practical inter'est for which considerable analysis can be made. Specifically, if the thrust acceleration is constant in magnitude and directed radially, tangentially, or circumferentially, then it is possible to obtain, at least partially, some quite interesting mathematical results.

A survey of a few more methods, which take into account exactly how you are planning to spiral out, that I find overall more accessible and useful in this case, is Petropoulos and Sims 2002, "A Review of Some Exact Solutions to the Planar Equations of Motion of a Thrusting Spacecraft". Note that they are most interested in the use of solar panels to power an ion thruster, allowing you to thrust in the most advantageous direction, which is very far from radial!

Answer 2

Partial answer only, I don' know what I'm doing exactly but I'll give this a try.

Update: Rockets work because their spectacular failures were scrutinized in the light of day rather than hidden. While this comment expresses concern that the math below might adversely confuse readers, I'm confident that to readers here it will be clear that it doesn't give the right answer, and that the technique of comparing a mathematical solution to a numerical solution serves as a good example for how an answer can be verified.

Math

Given initial acceleration at 1 AU is $8.17 \times 10^{-4} m/s^2$.

Tilt it at 45 degrees to make the thrust tangential, divide by $\sqrt{2}$ since its now oblique to the Sun, and account for drop off with distance from the Sun:

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

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

in the prograde direction (same direction as current velocity).

vis-viva equation:

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

Now

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

where $GM$ is the Sun's standard gravitational parameter $1.327 \times 10^{20} m^3/s^2$. The minus sign comes in because we know that contrary to first instinct, when we have an accelerating force in the prograde direction we counterintuitively decelerate b the same amount. This is cited in several other posts here as well, I'll look for some other answers to cite...

Rewrite and solve:

$$\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}$$

If we set $t(v_0) = 0$ in other words at time zero we are moving at orbital velocity $v_0 = \sqrt{GM/AU}$ we get 0.408 years.

Numerically

Interestingly when I try to simulate the same thing numerically, I get an escape time of 0.515 years! This surprises me for two reasons:

1. Velocity never reaches zero, meaning the math above is wrong!
2. And yet it's fairly close!

Conclusion

It will take roughly half a year to escape, but I don't have the full equation yet.

Further reading:

• Is it possible to reach the Sun without expending any fuel/reaction mass?
• What is the optimal angle for a solar-sail deorbit towards the Sun when radial thrust is included? (currently unanswered)

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()