Solar Sail Transfer Trajectory Design for Comet 29P/Schwassmann–Wachmann 1 Rendezvous

Abstract

The goal of this paper is to analyze the optimal transfer towards the periodic comet 29P/Schwassmann–Wachmann 1 of a solar sail-based spacecraft. This periodic and active comet is an interesting and still unexplored small body that has been regarded as an object of the Centaurs group. In this work, a classical (heliocentric) orbit-to-orbit transfer is studied from an optimal viewpoint, by finding the spacecraft trajectories that minimize the flight time for a given value of the solar sail characteristic acceleration, that is, the typical performance parameter of a photonic sail. In particular, the optimal Earth–comet transfer is studied both in a typical three-dimensional mission scenario and with a simplified two-dimensional approach, whose aim is to rapidly obtain an accurate estimation of the minimum flight time with a reduced computation cost. The numerical simulations illustrate the mission performance, in terms of the characteristics of the rapid transfer trajectory, as a function of the typical propulsive parameter and the solar sail thrust model.

1. Introduction

This paper investigates the performance of a solar sail-based spacecraft in a classical, orbit-to-orbit, three-dimensional heliocentric transfer from Earth to comet 29P/ Schwassmann–Wachmann 1 (SW1) [1,2,3]. Comet SW1, which was discovered on 15 November 1927 by German astronomers Arnold Schwassmann and Arno Wachmann, currently covers a nearly circular orbit between Jupiter and Saturn with a small inclination to the Ecliptic [4], as sketched in Figure 1. This periodic and active comet, whose orbital period is roughly 15 years, is an interesting and still unexplored small body that has been regarded as an object of the Centaurs group [5,6,7], that is, a celestial body originating from the Kuiper belt and then moved to the inner region of the Solar System. Thanks to its (relative) accessibility if compared with the Kuiper belt objects, comet SW1 has been proposed as a potential target for a robotic exploration mission [8]. In particular, the aim of the proposed missions was both to obtain in situ physical measurements of such a primordial celestial body [9,10] and to collect the particles ejected by the comet during its frequent outbursts [11,12].

In this paper, the capability of a spacecraft equipped with a photonic solar sail to perform a rendezvous of the SW1 orbit is investigated from an optimal viewpoint. In fact, the specific characteristics of the SW1 heliocentric orbit, which can be roughly approximated with a circular orbit of radius $6\mathrm{au}$, make the solar sail transfer a challenging problem, due to the fact that the maximum magnitude of the sail propulsive acceleration is inversely proportional to the square of the solar distance. Accordingly, and as also highlighted in the next section, just before the cometary rendezvous the propellantless propulsion system can give only a small fraction (about $1/6^2\simeq 3\%$) of the total thrust available just after the launch from the Earth, that is, when the Sun–spacecraft distance is roughly $1\mathrm{au}$. Indeed, a solar sail-based propulsion system is rarely considered in a heliocentric mission scenario which requires the thrust vectoring at a solar distance greater than $5\mathrm{au}$, that is, beyond the typical Jupiter distance from the Sun. However, the photonic solar sail can be still effective even beyond a radial distance of $5\mathrm{au}$, as highlighted by the numerical results discussed later in this paper.

In fact, the (rather) common assumption that a solar sail should be jettisoned, or the solar sail-induced propulsive acceleration should be neglected, when the solar distance is greater than $5\mathrm{au}$, originates from the interesting work by Carl Sauer [13] that analyzes the photonic solar sail performance in a Solar System escape (more precisely, the term used in [13] is ‘interstellar’) mission scenario. In particular, according to [13], in the last paragraph of p. 9 Sauer substantially affirms that “...the effectiveness of the sail in reducing flight time is minimal...” when the solar distance of the spacecraft is in the range $[5,10]\mathrm{au}$, so that the primary propulsion system (i.e., the solar sail) is jettisoned at $5\mathrm{au}$ “...to allow the acquisition of science data without possible interference from the sail”. This conclusion is, indeed, absolutely acceptable in a Solar System escape trajectory such as that analyzed in the second part of [13], because such a specific trajectory uses a pass close to the Sun (i.e., a sort of ‘solar photonic assist maneuver’ as named by Leipold and Wagner [14] and Dachwald [15] in their fine paper of 2004) to reach the required value of the hyperbolic excess speed. On the other hand, in a Solar System escape scenario the use of a solar sail-based propulsion system at a distance greater than $5\mathrm{au}$ is still possible, but the reduction in the total flight time (in that case, the flight times are of the order of decades) is marginal when compared to the increase in difficulty to actively change the sail attitude at such great solar distances.

Bearing in mind that a rendezvous with the active comet SW1 is therefore a challenging scenario for a (photonic) solar sail-based spacecraft, this work aims to illustrate the results of a preliminary, orbit-to-orbit, heliocentric trajectory design by considering the (simplified) thrust vector that models the performance of a flat sail of assigned characteristics. In this context, the paper takes advantage of both the analytical model and the ad hoc automatic routines already implemented by the authors to study this specific mission scenario. In this sense, the contribution of this work is not to propose new optimization techniques or to refine the mathematical model for the solar sail-based trajectory design, but to collect and comment on the results of a simulation campaign aimed at estimating the potentialities of such a propellantless propulsion system in a challenging mission scenario such as the rendezvous with the comet SW1.

In particular, taking the actual characteristics of the Earth and SW1 heliocentric orbits into account, the minimum flight time is evaluated by considering both an ideal force model (i.e., a model in which the sail is assumed as a flat, perfectly reflecting, mirror) and an optical force model (i.e., a model in which the thermo-optical characteristics of the reflecting film are taken into account) to describe the solar sail thrust vector during the deep space flight [16,17]. The numerical analysis presents a parametric study of the transfer time as a function of the sail typical performance parameter, i.e., the characteristic acceleration [18] defined as the maximum magnitude of the sail propulsive acceleration vector at a distance of 1 astronomical unit from the Sun.

Finally, the numerical results obtained through a simplified, but less computationally demanding, two-dimensional approach to the preliminary trajectory design are discussed in the last part of the work. In this context, bearing in mind the orbital characteristics (i.e., the value of the eccentricity and the inclination) of both the starting (Earth) and target (SW1) orbit, the spacecraft heliocentric transfer is approximated through a classical circle-to-circle orbit raising [19] and the mission performance is obtained by using the typical approach discussed in the recent literature [20]. The numerical simulations will show that, even with a non-negligible value of the target orbit inclination over the Ecliptic (the value of that orbital inclination is about 9 degrees), the two-dimensional circle-to-circle approximation of the orbit transfer gives both a flight time and a shape of the transfer trajectory which are very close to those obtained with a three-dimensional approach.

The paper is organized as follows. Section 2 briefly describes the mathematical model, that is, the three-dimensional spacecraft dynamics and the (rather classical) procedure which has been used to optimize the spacecraft trajectory during the interplanetary transfer. In particular, the aim of Section 2 is only to illustrates the main features of the mathematical model that, in its turn, has already been commented on in detail by the authors in the suggested references. Section 3 discusses the results of numerical simulations, in terms of a parametric study of the mission performance (i.e., the minimum flight time) as a function of the solar sail propulsive characteristics, and the results of a case study which refers to a specific value of the sail (reference) propulsive acceleration. A two-dimensional approach to the optimal trajectory design is then proposed in Section 4, whose aim is to obtain an approximate value of the optimal flight time with a reduced computation time. Finally, the concluding remarks are drawn in Section 5, as usual.

2. Mathematical Preliminaries and Models

The problem addressed in this section is to analyze the rapid transfer trajectory, in a three-dimensional heliocentric mission scenario, of a spacecraft whose primary propulsion system is a classical (photonic) solar sail of assigned characteristics. In this context, assume that the solar sail-based spacecraft initially (i.e., at time $t=t_0\triangleq 0$) covers a Keplerian (three-dimensional) orbit which coincides with the Earth’s orbit around the Sun, whose gravitational parameter is ${\mu }_{\odot }$. As is well known, such a specific scenario is consistent with a solar sail deployment just outside the Earth’s sphere of influence [21], when the spacecraft leaves the planet using a Keplerian parabolic escape trajectory.

In an orbit-to-orbit transfer, the spacecraft uses the propulsive acceleration given by the solar sail propulsion system to reach, at time $t=t_f>t_0$, a Keplerian orbit which coincides with the SW1 heliocentric trajectory. According to the common assumptions used in the preliminary phase of the trajectory design, during the transfer, that is, during the time interval $\Delta t\triangleq t_f-t_0\equiv t_f$, the spacecraft is subjected to the Sun’s gravitational acceleration and, possibly, solar sail-induced thrust. Moreover, in an orbit-to-orbit transfer the initial (or final) angular position of the spacecraft along the parking (or target) orbit is left free, so that the initial (or final) spacecraft true anomaly ${\nu }_0\triangleq \nu \left(t_0\right)$ (or ${\nu }_f\triangleq \nu \left(t_f\right)$) is obtained as an output of a trajectory optimization process that extremizes (i.e., optimizes) a given performance index. On the other hand, the remaining classical orbital parameters that define the parking (i.e., the Earth) and the target (i.e., the comet SW1) orbit are considered assigned, and they are reassumed in

Table 1. Classical orbital elements of the parking (Earth) and target (comet SW1) orbit used in the numerical simulations.

Body	a [au]	e	i [deg]	$\mathit{\omega }$ [deg]	$\mathsf{\Omega }$ [deg]
Earth	$1.001$	$1.594\times {10}^{-2}$	$3.022\times {10}^{-3}$	$302.978$	$159.864$
Comet SW1	$6.0499$	$4.4526\times {10}^{-2}$	$9.3619$	$51.4572$	$312.4009$

. In particular, this table gives the value of the semimajor axis a, eccentricity e, orbital inclination i, argument of perihelion $\omega$, and right ascension of the ascending node $\mathsf{\Omega }$, which are obtained through the well known JPL Horizons system.

Note that the comet SW1 orbital elements reassumed in Table 1 are evaluated at the beginning of August 2023. In fact, as discussed by Sarid et al. [22], the heliocentric orbit of comet SW1 is (strongly) perturbed by Jupiter’s gravitational attraction, and the simulations show that the Jupiter conjunction of year 2038 will (roughly) double the eccentricity of the comet orbit.

In a solar sail-based mission scenario, the time variation of the thrust vector components is typically calculated by minimizing the total flight time [23], that is, by studying the (rapid) transfer trajectory which maximizes the performance index J defined as

$$J\triangleq -\Delta t\equiv -t_f.$$

(1)

The procedure used in this work to maximize the performance index J and, then, to obtain the minimum flight time related to the rapid transfer trajectory, has already been validated in a number of different mission scenarios which involve the solar sail-based spacecraft transfer towards Solar System small bodies, as also illustrated in some recent authors’ works [24,25]. An interesting approach to the solar sail-based trajectory optimization, in missions toward Solar System small bodies, is also discussed in the works by Heiligers et al. [26] and Zeng et al. [27]. In particular, the procedure used in this work will be (briefly) illustrated at the end of this section, for the sake of completeness. The next subsection, instead, illustrates the mathematical model used to describe the solar sail-based spacecraft dynamics in the three-dimensional heliocentric scenario.

2.1. Spacecraft Dynamics

Paralleling the approach employed by the authors in the recent literature [25], the motion of the solar sail-based spacecraft during the heliocentric transfer is described using the modified equinoctial elements (MEOEs) introduced by Walker in 1985 [28,29]. In particular, MEOEs $\{p,f,g,h,k,L\}$ are elements valid for a generic osculating orbit, which are related to the classical orbital elements $\{a,e,i,\omega ,\mathsf{\Omega },\nu \}$ through the simple relationships [28,29]

$$\begin{array}{c}p=a\left(1-e^2\right),f=ecos(\omega +\mathsf{\Omega }),g=esin(\omega +\mathsf{\Omega }),h=tan(i/2)cos\mathsf{\Omega },\hfill \\ \hfill k=tan(i/2)sin\mathsf{\Omega },L=\nu +\mathsf{\Omega }+\omega \end{array}$$

(2)

while the inverse passage is obtained by using the equations

$$\begin{array}{cc}\hfill a& =\frac{p}{1-f^2-g^2}\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill e& =\sqrt{f^2+g^2}\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill i& =2arctan\sqrt{h^2+k^2}\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill sin\omega & =gh-fk,cos\omega =fh+gk\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill sin\mathsf{\Omega }& =k,cos\mathsf{\Omega }=h\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \nu & =L-\mathsf{\Omega }-\omega .\hfill \end{array}$$

(8)

Note that p is the semilatus rectum of the spacecraft osculating orbit, see the first of Equation (2), while the true longitude L is the fast varying variable which indicates the spacecraft angular position along the transfer orbit; see the last of Equation (2).

According to Betts [30], if one defines the spacecraft state vector $\mathit{x}$ as

$$\mathit{x}\triangleq {\left[p,f,g,h,k,L\right]}^{\mathrm{T}}$$

(9)

the spacecraft dynamics can be described by the following (compact) differential vectorial equation

$$\dot{\mathit{x}}=\mathbb{A}\mathit{a}+\mathit{d}$$

(10)

where $\mathit{a}\in {\mathbb{R}}^{3\times 1}$ is the solar sail-induced propulsive acceleration vector, while $\mathbb{A}\in {\mathbb{R}}^{6\times 3}$ and $\mathit{d}\in {\mathbb{R}}^{6\times 1}$ are defined as

$$\begin{array}{cc}\hfill \mathbb{A}& \triangleq \left[\begin{array}{ccc}0& \frac{2p}{1+fcosL+gsinL}\sqrt{\frac{p}{{\mu }_{\odot }}}& 0\\ sinL\sqrt{\frac{p}{{\mu }_{\odot }}}& \frac{\left(2+fcosL+gsinL\right)cosL+f}{1+fcosL+gsinL}\sqrt{\frac{p}{{\mu }_{\odot }}}& -\frac{g\left(hsinL-kcosL\right)}{1+fcosL+gsinL}\sqrt{\frac{p}{{\mu }_{\odot }}}\\ -cosL\sqrt{\frac{p}{{\mu }_{\odot }}}& \frac{\left(2+fcosL+gsinL\right)sinL+g}{1+fcosL+gsinL}\sqrt{\frac{p}{{\mu }_{\odot }}}& \frac{f\left(hsinL-kcosL\right)}{1+fcosL+gsinL}\sqrt{\frac{p}{{\mu }_{\odot }}}\\ 0& 0& \frac{\left(1+h^2+k^2\right)cosL}{2\left(1+fcosL+gsinL\right)}\sqrt{\frac{p}{{\mu }_{\odot }}}\\ 0& 0& \frac{\left(1+h^2+k^2\right)sinL}{2\left(1+fcosL+gsinL\right)}\sqrt{\frac{p}{{\mu }_{\odot }}}\\ 0& 0& \frac{hsinL-kcosL}{1+fcosL+gsinL}\sqrt{\frac{p}{{\mu }_{\odot }}}\end{array}\right]\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \mathit{d}& \triangleq \left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ \sqrt{{\mu }_{\odot }p}{\left(\frac{1+fcosL+gsinL}{p}\right)}^2\end{array}\right]\hfill \end{array}$$

(12)

The six scalar initial conditions associated with the vectorial differential Equation (10) model the spacecraft motion along the Earth’s heliocentric orbit at time $t_0$. Bearing in mind that in an orbit-to-orbit transfer the (initial) spacecraft angular position is left free, one obtains five initial conditions by combining Equation (2) and the data of the first row of Table 1. In that case, and according to [25], the result is given by the following equations

$$\begin{array}{cc}p\left(t_0\right)=p_0\triangleq 1.0005\mathrm{au},f\left(t_0\right)=f_0\triangleq -3.5430\times {10}^{-3}\hfill & ,g\left(t_0\right)=g_0\triangleq 1.5542\times {10}^{-2},\hfill \\ & h\left(t_0\right)=h_0\triangleq -2.4765\times {10}^{-5},k\left(t_0\right)=k_0\triangleq 9.0802\times {10}^{-6}\hfill \end{array}$$

(13)

The three components of vector $\mathit{a}$ which appear in Equation (10) coincide with the components of the propulsive acceleration vector in a typical radial–tangential–normal (RTN) reference frame of unit vectors $\{{\widehat{\mathit{i}}}_{\mathrm{R}},{\widehat{\mathit{i}}}_{\mathrm{T}},{\widehat{\mathit{i}}}_{\mathrm{N}}\}$. In particular, the RTN reference frame is sketched in Figure 2, in which $\mathit{r}$ is the Sun–spacecraft vector and $\widehat{\mathit{r}}\triangleq \mathit{r}/r$ is the radial unit vector, where r is the Sun–spacecraft distance. Note that ${\widehat{\mathit{i}}}_{\mathrm{R}}\equiv \widehat{\mathit{r}}$, the plane $({\widehat{\mathit{i}}}_{\mathrm{R}},{\widehat{\mathit{i}}}_{\mathrm{T}})$ coincides with the plane of the osculating orbit, and the normal unit vector ${\widehat{\mathit{i}}}_{\mathrm{N}}$ is aligned with the spacecraft specific angular momentum vector.

The analytical expression of the components of $\mathit{a}$ in the RTN reference frame depends on the solar sail thrust model which is adopted to describe the propulsive performance of such a propellantless thruster. In this context, assuming a flat sail, i.e., a solar sail with a rigid and flat reflecting membrane, a common choice in the preliminary trajectory design is to use the optical force model without degradation [31,32]. The interested reader should note that the solar sail thrust vector can be described through a more refined (and complex) mathematical model. In this sense, the work of Vulpetti [33], which can be considered one of the founders of the photonic solar sail astrodynamics, discusses the effect of wrinkles on the thrust vector characteristics. The optical force model, which is detailed in the classical McInnes’ textbook [34], allows the propulsive acceleration vector to be written, in a compact analytical form, as a function of the radial unit vector $\widehat{\mathit{r}}$ and the normal unit vector $\widehat{\mathit{n}}$, which is defined as the unit vector normal to the sail (nominal) plane in the direction opposite to the Sun [18,34]. Bearing in mind the approach described in [35], the compact expression of the propulsive acceleration vector is

$$\mathit{a}=a_c{\left(\frac{r_{\oplus }}{r}\right)}^2\left(\widehat{\mathit{r}}·\widehat{\mathit{n}}\right)\left[b_1\widehat{\mathit{r}}+\left(b_2\widehat{\mathit{r}}·\widehat{\mathit{n}}+b_3\right)\widehat{\mathit{n}}\right]$$

(14)

where $r_{\oplus }\triangleq 1\mathrm{au}$ is a reference distance, $a_c$ is the characteristic acceleration defined as the propulsive acceleration magnitude experienced by a solar sail whose reference plane is oriented perpendicular to the solar rays when $r=r_{\oplus }$, and $\{b_1,b_2,b_3\}$ are dimensionless (normalized) force coefficients where $b_1+b_2+b_3=1$, whose value depends on the optical characteristics of the reflecting membrane. In particular, assuming a membrane with a highly reflective aluminum-coated front side and a highly emissive chromium-coated back side [36], the values of the normalized force coefficients are reassumed in

Table 2. Normalized force coefficients for an ideal and an optical force model, assuming a (rigid and flat) sail membrane with a highly reflective aluminum-coated front side and a highly emissive chromium-coated back side.

Force Model	${\mathit{b}}_1$	${\mathit{b}}_2$	${\mathit{b}}_3$
Ideal	0	1	0
Optical	$0.0951$	$0.9109$	$-0.006$

, which also details the values of $\{b_1,b_2,b_3\}$ for an ideal force model [25].

Note that the term $a_c$ is the common solar sail performance parameter, whose value depends on the characteristics of both the solar sail and the whole spacecraft (including the value of the payload mass and other subsystems). For example, the upcoming NASA’s Solar Cruiser demonstration mission [37] has a characteristic acceleration of about $0.12\mathrm{mm}/{\mathrm{s}}^2$. Moreover, in Equation (14) the normal unit vector $\widehat{\mathit{n}}$ is the control term in the trajectory design, that is, bearing in mind the natural constraint $\parallel \widehat{\mathit{n}}\parallel =1$, one has two independent scalar control variables in the optimization problem, which is briefly described in the next subsection.

2.2. Trajectory Optimization

Fon an assigned value of the characteristic acceleration $a_c$, the optimal orbit-to-orbit transfer has been studied using an indirect approach [38] in which the Hamiltonian function $\mathcal{H}$ is obtained from Equation (10) as [25]

$$\mathcal{H}\triangleq \dot{\mathit{x}}·\lambda =\left(\mathbb{A}\mathit{a}\right)·\lambda +\mathit{d}·\lambda$$

(15)

where $\lambda \triangleq {\left[{\lambda }_p,{\lambda }_f,{\lambda }_g,{\lambda }_h,{\lambda }_k,{\lambda }_L\right]}^{\mathrm{T}}$ is the adjoint vector [39]. The time derivative of $\lambda$ is given by the classical Euler–Lagrange equations

$$\dot{\lambda }=-\frac{\partial \mathcal{H}}{\partial \mathit{x}}=-\frac{\partial (\mathbb{A}\mathit{a}·\lambda )}{\partial \mathit{x}}-\frac{\partial (\mathit{d}·\lambda )}{\partial \mathit{x}}$$

(16)

which gives six first-order, non-linear, differential equation in the variables $\{{\lambda }_p,{\lambda }_f,$${\lambda }_g,{\lambda }_h,{\lambda }_k,{\lambda }_L\}$, whose explicit expressions are omitted here for the sake of brevity.

During the transfer, the optimal control vector, that is, the optimal sail attitude in terms of time variation of $\widehat{\mathit{n}}$, is obtained by using Pontryagin’s maximum principle, that is, by finding the value of the normal unit vector that maximizes (at any time) the Hamiltonian function given by Equation (15). In this respect, a numerical approach is used to calculate the two (independent) scalar variables that maximize the value of $\mathcal{H}$ and the procedure is illustrated in the recent authors’ work [25].

The equations of motion (10) and the Euler–Lagrange Equation (16) form the two-point boundary problem (TPBVP) which is completed by 13 boundary conditions, because the (optimal) flight time $t_f$ is an output of the optimization process. In particular, five (initial) conditions are given by Equation (13), while five (final) conditions are obtained by enforcing the rendezvous with comet SW1 orbit. In this respect, using the data from Table 1 and Equation (2), one obtains

$$\begin{array}{c}p\left(t_f\right)=p_f\triangleq 6.0379\mathrm{au},f\left(t_f\right)=f_f\triangleq 4.4425\times {10}^{-2},g\left(t_f\right)=g_f\triangleq 2.996\times {10}^{-3},\hfill \\ \hfill h\left(t_f\right)=h_f\triangleq 5.5214\times {10}^{-2},k\left(t_f\right)=k_f\triangleq -6.0465\times {10}^{-2}\end{array}$$

(17)

Finally, the last three scalar constraints required to complete the TPBVP are obtained from the transversality condition [39] as

$${\lambda }_L\left(t_0\right)={\lambda }_{L_0}\triangleq 0,{\lambda }_L\left(t_f\right)={\lambda }_{L_f}\triangleq 0,\mathcal{H}\left(t_f\right)={\mathcal{H}}_f\triangleq 1$$

(18)

The TPBVP has been solved using an in-house procedure based on the multiple shooting method and the use of genetic algorithms. In particular, the differential equations have been integrated in double precision using a variable order Adams–Bashforth–Moulton solver [40] with absolute and relative errors of ${10}^{-10}$. The results of this numerical optimization process are illustrated in the next section, which contains both a parametric analysis of the orbit-to-orbit transfer as a function of the value of the characteristic acceleration and a detailed case study for a specific value of $a_c$.

3. Numerical Simulations

The optimization procedure discussed in the previous section has been applied to evaluate the rapid Earth–comet SW1 transfer as a function of the spacecraft characteristic acceleration value in the range $a_c\in [0.5,2]\mathrm{mm}/{\mathrm{s}}^2$, with a step of $0.1\mathrm{mm}/{\mathrm{s}}^2$; see the red circles in Figure 3 for example.

Bearing in mind the current solar sail technology level [16,23], the range below the (usually considered as canonical [34]) value of $1\mathrm{mm}/{\mathrm{s}}^2$ is consistent with a medium–high performance solar sail, which hopefully will be available in the near future. On the other hand, the upper part of the selected range of $a_c$ models the behavior of a (more futuristic) high performance sail [41,42,43,44], which uses advanced materials [45,46] to improve the propulsive performance. In order to evaluate the impact on the transfer performance of the sail thrust characterization, the orbit-to-orbit transfer has been optimized by considering both the ideal and the optical force model, using the normalized coefficients $\{b_1,b_2,b_3\}$ reassumed in Table 2. The results, in terms of optimal flight time $\Delta t$ as a function of $a_c$, are given by Figure 3 for the two force models.

Figure 3 clearly shows that the (optimal) flight time of a medium–high performance solar sail is well above 10 years, with a value of $\Delta t\simeq 25\mathrm{years}$ when $a_c=0.5\mathrm{mm}/{\mathrm{s}}^2$ for both the force models considered in the numerical simulations. On the other hand, an orbit-to-orbit transfer time below 10 years can be reached only through a high performance solar sail. In this context, for example, a sail with $a_c>1.8\mathrm{mm}/{\mathrm{s}}^2$ is required when an optical force model is used to describe the thrust vector. This aspect is due to the combination of a high value of the semimajor axis and a low eccentricity of the comet SW1 heliocentric orbit, which gives a perihelion (or aphelion) distance of about $5.78\mathrm{au}$ (or $6.32\mathrm{au}$). In fact, a classical solar sail is quite ineffective at such distances from the Sun, because the maximum magnitude of the solar sail propulsive acceleration $max\parallel \mathit{a}\parallel$ varies as $1/r^2$, as described by Equation (14). In this specific mission scenario one obtains that the dimensionless ratio $max\parallel \mathit{a}\parallel /a_c$ is in the range $[0.025,0.03]$ when the spacecraft completes the rendezvous, as illustrated by the curve of Figure 4, which describes the variation of the ratio $max\parallel \mathit{a}\parallel /a_c$ as a function of the spacecraft (final) true anomaly ${\nu }_f$ along the comet orbit. In other terms, at the end of the mission the solar sail-induced propulsive acceleration magnitude is only a small fraction (i.e., a value below $3\%$) of the initial maximum value, which roughly coincides with $a_c$ by construction.

The numerical simulations show that the spacecraft with a medium–high performance solar sail completes several revolutions during the optimal orbit-to-orbit transfer. This point is highlighted by Figure 5, which shows the variation in true longitude $\Delta L$ during the transfer. Note that the integer part of the ratio $\Delta L/\left(2\pi \mathrm{rad}\right)$ indicates the revolutions around the Sun completed by the spacecraft before reaching the comet SW1 heliocentric orbit. The value of $\Delta L$ is consistent with the curves reassumed in Figure 6, which shows the ecliptic projection of the optimal transfer trajectory for some values of $a_c$, using an optical force model.

Figure 3 also shows, as expected, that the performance of the ideal force model (in terms of flight time) is better than that obtained through an optical force model. However, the relative difference between the (minimum) flight time $\Delta t$ obtained through the two force models is below $10\%$, as highlighted in Figure 7.

The numerical solution of the TPBVP gives, as secondary outputs, the (optimal) value of the true anomaly both at the beginning (${\nu }_0$) and the end (${\nu }_f$) of the transfer, as sketched in Figure 8. Note that the value of the pair $\{{\nu }_0,{\nu }_f\}$ is useful to find a launch window for the rapid transfer, as a function of the value of the characteristic acceleration.

Case Study

This subsection gives a detailed calculation of the characteristics of the optimal transfer trajectory for specific value of $a_c$. In this context, the (canonical) value of $a_c=1\mathrm{mm}/{\mathrm{s}}^2$ has been selected to describe the time variation of both the spacecraft osculating orbital parameters and the sail attitude during the orbit-to-orbit transfer. In this case, the optimal flight time is $\Delta t=5306\mathrm{days}$, that is, about $14.5\mathrm{years}$; see also Figure 3b. The time variation of the classical orbital parameters of the spacecraft osculating orbit is shown in Figure 9, where one can observe that the change in orbital inclination is substantially performed in the first part of the transfer, that is, when the Sun–spacecraft distance is still in the order of one astronomical unit. The sketch of the optimal transfer is reported in Figure 10.

Finally, the optimal control law, in terms of the time variation of the component of the normal unit vector $\widehat{\mathit{n}}$ in the RTN reference frame, is reported in Figure 11. Note that the three components of $\widehat{\mathit{n}}$ are roughly constant for $t\ge 5\mathrm{years}$, that is, the sail attitude (with respect to the RTN reference frame) is nearly fixed for about $60\%$ of the transfer. In particular, such an attitude maximizes the transverse component (i.e., the component perpendicular to the radial direction) of the propulsive acceleration $\mathit{a}$. This is an interesting result, which is consistent with the model described in the next section, where a two-dimensional approach is used to obtain an accurate approximation of the optimal transfer performance with a reduced computation time.

4. Two-Dimensional Simplified Approach

The (very) small value of the eccentricity of both the Earth and comet SW1 heliocentric orbit, together with the small relative orbital inclination (slightly above $9\mathrm{deg}$, see data in Table 1), suggest the use of a simplified mathematical model to study the optimal orbit-to-orbit transfer. In fact, an Earth–comet SW1 mission is the typical scenario that can be approximated through a classical circle-to-circle orbit raising, in which the celestial bodies’ (heliocentric) orbits are assumed as circular and coplanar. In that case, the spacecraft dynamics can be described with a (simplified) mathematical model which uses the classical spacecraft equations of motion in the polar (heliocentric) reference frame sketched in Figure 12. In this respect, the dynamics model and the optimization procedure are thoroughly described (and commented on) in a number of the authors’ works [19,47], and they are omitted here for the sake of brevity.

In particular, the radius of the initial (or final) circular orbit is assumed as $r_0=r_{\oplus }=1\mathrm{au}$ (or $r_f=6.0499\mathrm{au}$, i.e., a value that coincides with the semimajor axis of the comet SW1 orbit), while the initial polar angle is set to $\theta \left(t_0\right)=0$ for symmetry reasons; see Figure 12. Using such a two-dimensional approach to trajectory optimization, and considering the range of characteristic acceleration $a_c\in [0.5,2]\mathrm{mm}/{\mathrm{s}}^2$, the minimum flight time for an ideal or optical force model has been reassumed in Figure 13 as a function of $a_c$. A comparison between the curves of Figure 13 and that of Figure 3 (which refer to a complete three-dimensional scenario), indicates that the results obtained through a two-dimensional approach approximate very well the values obtained with a more realistic model. This aspect is confirmed by Figure 14, which shows that the percentage difference between the minimum flight time obtained through a three-dimensional (${\Delta t|}_{3\mathrm{D}}$) and a two-dimensional (${\Delta t|}_{2\mathrm{D}}$) approach is below $4.6\%$ for both the two force models.

The solution of the optimal control problem also gives the value of the final polar angle $\theta \left(t_f\right)={\theta }_f$, which can be used to evaluate the number of complete revolutions around the Sun during the transfer. In this case, the numerical results are reassumed in Figure 15, while Figure 16 shows a set of optimal transfer trajectories (in polar form), when an optical force model is considered in the simulations. Bearing in mind that $\theta \left(t_0\right)=0$, one observes that the trajectories of Figure 16 are consistent (and in some cases substantially coincident) with the trajectories reported in Figure 6. This is an important result, because it confirms a posteriori the usefulness of the two-dimensional approach, which requires a computation time of about two orders of magnitude less than that related to the three-dimensional model.

5. Conclusions

A solar sail-based rendezvous with comet 29P/Schwassmann–Wachmann 1 requires a flight time of about two decades when a medium–high performance sail is considered in the orbital simulations. Accordingly, such a complex mission scenario can be considered a typical application of a high performance sail, which could give a flight time less than ten years, even with a realistic optical force model. In this respect, the numerical results discussed in this paper can be used to obtain a fast estimation of the characteristics of the optimal transfer trajectory with a reduced computation time. A possible extension of this work is to consider a cometary flyby mission, instead of a rendezvous scenario, in order to reduce the required flight time for a given value of the characteristic acceleration. In fact, a flyby mission has a less constrained terminal condition, so that the flight time is surely lesser than the time obtained in the paper. Accordingly, for example, if the scientific mission requires a cometary flyby and not an in situ study, the flight time can be significantly less than the decades obtained in the previous section, the value of the characteristic acceleration being the same.