Uppsala University

(1)

Uppsala University

Student Project

Supervisor: Dr. David Andrews

Full Physics Orbital Simulation around Comet Didymos for CubeSat Mission APEX on HERA, ESA

Author:

Jonas Zbinden

(2)

Abstract

We investigate the dynamic environment of the double asteroid system (65803) Didymos with its satellite Didymos B for stable trajectories of more than 30 days for a model cubesat called APEX, foreseen to visit Didymos with the HERA mission, to be launched in 2023. We find semi-stable orbits with lifetimes up to 15 days under the influence of solar radiation pressure. In comparison, in the absence of solar radiation pressure, orbit trajectories reach lifetimes of up to 90 days. We present the distribution of orbital stability for positions of initial deployment and the related orbital lifetime for inclinations of i = 0◦, 6◦ and 30◦ to the ecliptic plane of the sun, to investigate expected seasonality of the Didymos system. The influence on trajectories through solar radiation pressure on the arrangement of the solar panels of APEX is compared to trajectories without solar radiation pressure. It becomes clear that the influence of solar radiation pressure is the main variable for stability of trajectories and must be further investigated with great care. We also discuss briefly the operation of APEX in relation to orbital corrections and the limitations for such corrections.

1 Introduction

(3)

surfaces like solar panels. We adopt their approach to the cubesat mission APEX, designed by the APEX consortium and compare their results to our results in which we use the ac-tual shape, surface properties and mass of the cubesat. We will also discuss how one could plan course corrections and therefore navigate safely in even closer proximity of Didymain or Didymoon, and possibly "land" on Didymoon.

2 Theory

2.1 Equation of motion

The equation of motion of a cubesat (APEX) in the restricted three body problem in an intertial frame of reference with ~r denoting the position of APEX relative to Didymain is given by: ¨ ~ r = ¨~rSRP(~r, t) − GMDidM r3 ~r + ¨~rdeg2(~r, t) − GMDidm r3_Didm rDidm~ (1) where on the left hand side we have the acceleration ¨~r on APEX and on the right hand side we have the different force contributions acting on APEX, normalized to the mass of APEX. We have in the following order, the solar radiation pressure (SRP), the leading order term of the gravitational force of Didymain, denoted with DidM , the secondary or higher order terms of the gravitational force of Didymain, that arise due to its oblate symmetric ellipsoidal shape and the leading order term of the gravitational force of Didymoon, denoted with Didm and rDidm~ = ~r − ~rorb with ~rorb the vector from Didymain to Didymoon. The

mass of APEX solely shows up in the solar radiation pressure term, no other dissipative forces are considered.

2.2 Solar radiation pressure

The solar radiation pressure for each surface (sur) can be represented with ¨ ~rsur = −P(1AU d )2 A ms/c cos (θ)[(1 − )~e+ 2 cos (θ)~n] (2)

where P = 4.56 · 10−6 N/m2 is the solar radiation pressure at 1 AU , d is the distance

from the surface of the sun, A is the surface area, m_s/c is the mass of the spacecraft, θ is illumination angle on to the surface, which means the angle between the incoming radiation vector ~e and the normal vector of the surface ~n. is a material constant describing the

(4)

Figure 1: Dimensions and orientation of APEX cubesat system.

(5)

of Didymoon introduces a seasonal variability in solar radiation pressure. Since there is no fixed schedule available yet for the mission of HERA, we have not quantified these seasonal effects, but rather investigated the influence of different illumination angles with respect of the orbital plane of Didymoon and fixed distance from the sun at 1 AU. In case of APEX trajectory leading behind Didymain or Didymoon, we included appropriate shadow functions to turn off the solar radiation pressure, as soon as the position of APEX fulfils the shadow criteria. Obviously there is the potential of erroneous solar radiation influence, since between two integration steps APEX could travel significantly far into a shadow region. However, we assume the integration stepsize normally to be smaller than the precision at which we simulate the shadow for which we assume spherical bodies. All these simplifications would need to be considered again once APEX has been released at Didymos and we had a better picture of the shape and motion of Didymos. The shadow criterias were deduced from [5], leaving out the treatment of the Penumbra region due to the small angle of the penumbra cone f1:

sin (f1) = (R+ RDidM)/s≈ R/s≈ 10−8 (3)

⇒ f1 ≈ 10−8. (4)

with s = d the distance from the sun.

(6)

2.3 Gravitational Potential of Didymos

Since gravitation is a conservative force we can express forces from Didymos acting on APEX as the gradient of the gravitational potential of Didymain and Didymoon. For both bodies we have written down the leading order term of the gravitational potential in the equation of motion (1). On top of that we include the higher order terms of the expansion of the graviational potential of Didymos up to second order, following the procedure of Damme et al. [4]. They investigated the contribution of higher order terms and find that the effect of higher than second order was not noticeable. The full expansion of the gravitational potential in spherical harmonics is represented by [5]

U (r, φ, λ) = GM r ∞ X n=1 n X m=1 Rn

rnPnm(sin (φ))(Cnmcos (mλ) + Snmsin (mλ)). (5)

λ = π₂ − θ where θ, φ are spherical coordinates. Only the terms with n, m ∈ {0, 2} yield non-zero contributions for an oblate symmetrical ellipsoid. The leading order term n = 0 has already been excluded here. Snm = 0 identically, due to the symmetry of Didymain.

R is the radius Pnm are the associated Legendre polynomials. The unnormalized spherical

harmonic coefficients were again adopted from Damme et al. [4]. We are left with U2(r, φ, λ) = − 1 r3GMDidMR 2 1 2C2,0 3 sin (λ)2− 1 + 3C2,2cos (λ)2cos (2φ). (6)

The minus in front of U₂ is a chosen by convention, since we use the convention for the force ~

F = −∇U . The complete gravitational potential of Didymos can be transformed into the co-rotating frame of Didymoon around Didymain. The angular velocity can be calculated with the orbital period which we know from Damme et al. [4]: ω = 2π_T , T = 11.92 h. The following coordinates refer now to coordinates in the corotating frame:

Ucorot(r, φ, λ) = U − 1 2ω 2_r2 ₍₇₎ = −GMDidM r − GMDidm rDidm + U2− 1 2ω 2_r2 ₍₈₎

where rDidM = r and rDidm= |~r−~rorb|. The gravitational potential of Didymos is presented

(7)

L1 L2 L3

L4?

L5?

Figure 3: The gravitational potential clearly reveals the positions of L1, L2 and L3, L4 and L5 seam not resolvable, however, as shown later there exist stable orbits from L4, L5 and L3.

(8)

Asteroid Didymos (main) Didymos (moon) Mass [kg] 5.24 ×1011 3.45 ×109 Mean Radius [km] 0.385 0.075 Axes a, b, c [km] 0.395, 0.39, 0.37 -C2,0 unnormalized -0.023 -C2, unnormalized -0.0013 -Spin period [h] 2.26 -Mutual orbit Orbital period [h] 11.92 Semi-major axis [km] 1.178 Eccentricity 0.05 Orbital pole A (300, -60) [◦] Orbital pole B (310, -84) [◦]

Table 1: Data used for the Didymos system in our orbit calculations. For original reference refer to table 1 in Damme et al. [4], where the data presented here is reprinted from. To get the force we computed the gradient of U_corot respectively of U while conversion to the co-rotating frame is not necessary anymore, since we solve the differential equation in an inertial frame. The gradient for U2 was computed by hand and checked with Mathematica.

2.4 Differential equation solver and parallel computation

To solve the differential equation we used thescipy.integrator.odesolver. We used the inte-grator "dop853" which uses an excplicit Runge-Kutta order 8(5,3) scheme with automatic stepsize control following the "dopri" method, developed by Dormand and Prince. It uses an error estimate from the difference between the fourth and fifth order solution to adjust its next stepsize. [6] To reduce computation time we implemented an automatic process distributor for parallel computation, called joblib [7].

3 Results

First investigations were performed to find stable orbits around Didymos. To assess which orbits can be regarded as stable or unstable, the following abort criteria were used:

r < RDidM (9)

rDidm < RDidm and (10)

r > 1000 × a (11)

(9)

positions coarse, the trajectories were computed for up to 90 days. The initial coordinates in relation to the trajectorie’s lifetime are presented below.

−2000 −1500 −1000 −500 0 500 1000 1500 2000 x-coordinate [m] 2000 1500 1000 500 0 500 1000 1500 2000 y-c oo rdi na te [m ]

Orbit Lifetime from initial conditions of APEX in corotating frame

0 5 10 15 20 25 30 Or bit Li fet im e

Figure 4: Orbital lifetimes for APEX with different deployment positions depicted as the coloured dots. The color indicates the orbital lifetime of an orbit from that particular position. The initial velocity is computed for each point according to its position and matches the co-rotation of Didymoon around Didymos with influence from solar radiation pressure.

(10)

Orbit Lifetime from initial conditions of APEX in corotating frame

0 5 10 15 20 25 30 Or bit Li fet im e

Figure 5: For comparison, in case of no influence from solar radiation pressure, the orbital lifetimes were generally longer and there we found more stable orbits with orbital lifetimes of at least 30 days. Possible orbits of up to 90 days stability were observed as well.

(11)

−2000 −1500 −1000 −500 0 500 1000 1500 2000 x-coordinate [m] −2000 −1500 −1000 −500 0 500 1000 1500 2000 y-c oo rdi na te [m ]

Orbit of APEX in corotating frame Circ lar Orbit Didymoon

Solar radiation press re on

(a) Shape of orbit APEX released from L4

Solar radiation press re on

(b) Shape of orbit APEX released from L5

Figure 6: Comparison between the shape of orbits, both released at Lagrange points L4 and L5. It can not be distinguished from which Lagrange point APEX was released and show how the space craft travels between the two Lagrange points semi-stable, eventually colliding with the main body.

Solar radiation press re off

(a) Shape of orbit APEX released close to L4

Solar radiation press re off

(b) Shape of orbit APEX released from L5

(12)

The stable orbits are concentrated inside or outside the binary system. Both prograde and retrograde orbits in comparison to the rotation of Didymos could be found. Additionally we found more chaotic orbits following the distribution of L3, L4 and L5. This stands in agreement with the findings of Damme et al. [4], except of the somewhat more enlarged area the orbits cover. This can be explained by the two following reasons. On one hand our shape and mass of APEX differ from the assumptions of Damme et al. [4]. Additionally we observed that for the same initial conditions, the orbits do not necessarily lead to the same result. This might be explained by numerical dissipation like stepsize, internal errors like storage of variables or by how precise the ODE solver works. This has nothing to do with chaos, however, on top of our observations in comparison to Damme et al. [4], we observe highly chaotic orbits if the influence of solar radiation pressure is included.

To investigate the influence of different inclinations of the orbital plane of Didymos towards the ecliptic of the sun, we used the orbital poles adopted from table1in Damme et al. [4] and included the inclinations 30◦and 6◦. The results are presented in the following figures.

Orbit Lifetime from initial conditions of APEX in corotating frame

0 5 10 15 20 25 30 Or bit Li fet im e

(13)

Orbit Lifetime from initial conditions of APEX in corotating frame

0 5 10 15 20 25 30 Or bit Li fet im e

Figure 9: Orbit lifetimes with 30 ◦ inclination. The number of orbits lasting longer than 10 days and the maximum orbit lifetime are similar to no inclination to the ecliptic of the sun. The number of orbits lasting longer than 10 days even doubled.

The inclination has an influence on the orbital lifetime and for greater inclination even a stabilizing one. Qualitatively the regions of safer deployment for APEX do not change. Additionally, in case of i = 30◦ we found one more stable orbit closely above Didymoon.

(14)

4 Discussion and Conclusion

We have shown the existence of some areas with orbital lifetimes up to 15 days. Although we have integrated the e. o. m. for up to 90 days, no trajectories were found that could survive longer than 17 days under the influence of solar radiation pressure. Without solar radiation pressure, trajectories respectively. orbits with orbital lifetime of > 90 days could be found. The influence of radiation pressure is the main perturbation for our system. Looking at the shape of APEX and having in mind its change of orientation along an orbit with respect to the sun, the perturbation of the trajectory of APEX is not predictable for a general setup, but must be evaluated on sight for the true orientation of APEX. A solution to "stabilize" the influence of solar radiation pressure could be turning APEX in a "stand-by" position, always facing the same side towards the sun. Thereby the force due to the illumination on the solar panels of APEX could be minimized. In any case, this issue needs further investigation.

For all our trajectory calculations we used a time step size of dt = 100 s. The ODE solver adjusts the stepsize automatically to comply with the conversion criteria. Nevertheless, we saw minor changes in precision for lower dt. Another issue is that we do not see more stable orbits in contrast to the work of Damme et al. [4]. This can be explained by the different shape of APEX, which might lead to chaotic behaviour under arbitrary orientations of APEX compared to the sun. On the other hand it can of course also be the case that we did not find all regions of stable orbits due to a too coarse grid of initial conditions or maybe the choice of dt. However the convergence criteria used in the ODE solver should normally solve the issue of choosing dt too big. We also observed a sensitivity in the precision of the calculation, since runs with the same initial conditions could lead sometimes to different results. If the precision fails on the computation side (storage of variables) or the ODE solver side, could not be answered and requires further investigation.

One issue we have not looked at, is if we might have missed collisions due to the choice of dt. This is discussed as well by Damme et al. [4] and they include a polygon chain for each position of their model satellite to check if they missed a collision. In a future investigation this should be included in our tool as well.

(15)

can say that a bad choice of initial position can not be corrected by a velocity correction of up to 1 ms−1 in most cases. Interestingly during the velocity correction analysis we found fluctuations in v_z, although we assumed no inclination towards the ecliptic of the sun and did not apply velocity corrections in z. This can only be explained by technical artefacts from the computation during the integration of our solution.

We conclude that there is more work to do to investigate the different influences on the simulation results, both from physical sources like the orientation of Didymos to the sun and the distance from the sun, the orientation of APEX compared to the sun as well as numerical influences like the sensitivity to stepsize, and numerical accuracy. However, to address these issues more time and computation power is needed than intended or available for this project.

(16)

conservative forces are time independent. The solar radiation pressure would be regarded as flipped in sign, since it is a non conservative force due to its partially inelastic nature, namely, absorption of photons. This is physically not completely correct since flipping the arrow of time in case of inelastic collisions results in photons jumping out of the material and taking energy from it. However for our purposes, the assumptions made previously would be enough. Only the initial velocity vector would flip in sign ~v ⇒ −~v which is similar to "throwing" APEX from the surface of Didymoon and observing in which orbit around the Didymos system it ends up in. This immediately implies that for a chosen landing site, APEX must be able to end up in such an orbit. Due to its rather chaotic behaviour, this might be as difficult to achieve as to simply calculate all possible orbital corrections from a certain position which would lead to a controlled crash-landing on Didymoon.

Our software tool provides many capabilities to analyze the dynamic environment of Didy-mos for a cube sat like APEX with a high cost in computational capacity.

References

[1] A.F. Cheng, P. Michel, M. Jutzi, A.S. Rivkin, A. Stickle, O. Barnouin, C. Ernst, J. Atchison, P. Pravec, and D.C. Richardson. Asteroid impact and deflection assessment mission: Kinetic impactor. Planetary and Space Science, 121:27 – 35, 2016.

[2] Patrick Michel, A. Cheng, M. Küppers, P. Pravec, J. Blum, M. Delbo, S.F. Green, P. Rosenblatt, K. Tsiganis, J.B. Vincent, J. Biele, V. Ciarletti, A. Hérique, S. Ulamec, I. Carnelli, A. Galvez, L. Benner, S.P. Naidu, O.S. Barnouin, D.C. Richardson, A. Rivkin, P. Scheirich, N. Moskovitz, A. Thirouin, S.R. Schwartz, A. [Campo Bagatin], and Y. Yu. Science case for the asteroid impact mission (aim): A component of the asteroid impact and deflection assessment (aida) mission. Advances in Space Research, 57(12):2529 – 2547, 2016.

[3] Hera mission & cm16 lessons learned. https://www.cosmos.esa.int/ documents/336356/1503750/SMPAG09_HERA_Carnelli_2017-10-11.pdf/

c2c0e63d-256f-6335-1229-7e1fea4cd733. Accessed: 2020-06-08.

[4] Friedrich Damme, Hauke Hussmann, and Jürgen Oberst. Spacecraft orbit lifetime within two binary near-earth asteroid systems. Planetary and Space Science, 146:1 – 9, 2017. [5] O. Montenbruck and E. Gill. Satellite Orbits: Models, Methods, and Applications.

Physics and astronomy online library. Springer Berlin Heidelberg, 2000.

[6] Syvert P. Nørsett Ernst Hairer, Gerhard Wanner. Solving Ordinary Differential Equa-tions I, volume 2. Springer, Berlin, Heidelberg, 1993.