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

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)
• I an new to Stack Exchange @peter nazarenko how did you get this great re-formatting. Is there a manual for how to do math symbols? – John Davies Dec 10 '20 at 15:57
• Yes, there is a little help for format, including math symbols. You can press the "(?)" mark at the upper right corner of the editing subwindow, and then press "Advanced help" hyperlink and then scroll down (math formatting described near the bottom of the Advanced Help page, titled as "LaTeX"). Welcome to Space Exploration StackExchange and Good Luck! – Peter Nazarenko Dec 10 '20 at 16:48
• In order to escape the sun, the sail must be pointed not directly outward, but at an angle, to provide a forward acceleration. If the sail is pointed directly outward, the effect on the spacecraft's orbit would be the same as if the Sun's mass decreased to $(1-\frac{8.17\times 10^{-4}}{5.9\times 10^{-3}}) = 0.862$ of its value. The orbit will simply become elliptical. – Litho Dec 11 '20 at 10:59
• @Litho: you're a genius. I was busily writing a comment saying it would need to be attacked numerically and ... no, obviously not. – tfb Dec 11 '20 at 11:02
• @litho thank you, your comment makes perfect sense. It would be an interesting problem to work out what the optimum angle of the sail would be for the spacecraft to get greatest distance from the Sun in the shortest time. – John Davies Dec 11 '20 at 11:54

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!

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

$$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.

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)

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

• – uhoh Dec 12 '20 at 3:42
• I think this answer should be removed until the issues with the calculation are solved. In the current form it might be confusing to readers. – asdfex Dec 13 '20 at 9:43
• @asdfex I think that the introduction "Partial answer only, I don' know what I'm doing exactly..." is sufficient warning to avoid that, and SE works on the principle that answers without up votes are less likely to be correct. This post explains that there is some dissonance between the first attempt at an analytical and the numerical reality, and serves as a warning that an analytical solution will have to do better. Science and math build from failed attempts. Rockets work because their spectacular failures were scrutinized in the light of day rather than hidden. – uhoh Dec 13 '20 at 22:19
• @asdfex but I'll fortify the first sentence accordingly. – uhoh Dec 13 '20 at 22:21