Rendezvous of a continuous low-thrust cubesat with a satellite

We give in this paper a new method optimizing continuously the direction of a low continuous thrust to have rendezvous of a cubesat with a given satellite. We study the rendezvous of a cubesat B, with a satellite A, on an elliptic orbit around the Earth. The cubesat B, with 3 kg mass and low-thrust motor, is on a circular orbit around the Earth. The rendezvous is obtained by making the coincidence between the angular momentum vector and the eccentricity vector of B respectively with those of A, then by adjusting the position of B, at the initial instant of thrusting or by convenient choice of this instant. The thrust acceleration is supposed to be constant. The thrust is continuous and its direction is continuously oriented to optimize the time necessary to make the coincidence between the two orbits. A numerical application is given where A is the satellite ‘Almasat-1’.


Introduction
The low-thrust spacecraft is very useful for the exploration of the space around the Earth, or around any planet in our solar system. For a given mission, the rendezvous, by continuous low-thrust, is obtained after a long time, but it needs a small quantity of energy. The low-thrust can be produced by a motor which use electrical ionisation. We suppose that the direction of the thrust is controlled at any instant.
Many researches 5 use low-thrust to perform rendezvous with asteroids or to make transfers to the Moon or around Mars, or to have optimal Earth-Moon trajectories 6 . More general information about low-thrust studies can be obtained from [14][15][16][17][18][19]. In these studies, the thrust is not continuous, but in this paper we use a continuous low-thrust which is continuously oriented to optimize the time necessary to obtain a rendezvous of a cubesat with a given satellite.
We consider a small low-thrust spacecraft, a cubesat B having a mass of 3 kg on a given elliptical orbit around the Earth. The cubesat B is the interceptor. The target is a given non manoeuvrable satellite A having a known position at an instant t 1 on an elliptical orbit around the Earth. The motor of B produces a constant continuous thrust whose direction is not constrained. We switch on B's engine at a given instant t 1 in order to make the rendezvous with A, by varying adequately the thrust direction. This amounts to coinciding, at a time t 2 >t 1 , the position and the velocity of B with those of A. We will see how the satellite B can intercept the satellite A. In fact the orbit of the interceptor B is Keplerian when the thrust is stopped, but when the thrust is continuous it is never Keplerian; it is an osculating orbit which is constantly deformed (like a work of potter slowly shaping an amphora) in order to coincide with the orbit of A at the instant t 2 .
We propose a numerical method to achieve the rendezvous of B with A.

Equations and notations
We use an inertial direct frame Oxyz whose origin is the centre O of the planet, the Earth for example. An Oz axis passes through the North Pole. An Ox axis in the plane of the equator, points to a fixed inertial direction of space, for example the vernal equinox. The classical parameters defining the orbit of a satellite in this inertial reference frame are: a Semi-major axis of the ellipse, e Eccentricity of the orbit, i Inclination of the orbit relative to a reference plane, for example the equator.
Ω Longitude of the ascending node, angle measured in the plane of the equator from Ox axis to the line of intersection with the plane of the orbit.
ω Argument of pericentre, angle measured in the plane of the orbit from the ascending node to the pericentre.
M Mean anomaly, angle measured from the pericentre to the radius vector of a fictitious satellite in a circular orbit of same major axis. This angle helps calculating the true anomaly giving the instantaneous position of the satellite on its orbit.
In the first step, we obtain the coincidence of the B orbit plane with the A orbit plane, by adjusting the angular momentum vector h B  of B with the angular momentum vector h A  of A. That is equivalent to make Ω B and i B coincide with Ω A and i A respectively.
In the second step, we obtain the coincidence of the B orbit with the A orbit, by adjusting the eccentricity vector e B  of B with the eccentricity vector e A  of A. That is equivalent to make e B and ω B coincide with e A and ω A respectively.
We assume that the thrust acceleration, or thrust per unit mass, T  (t)/m (t), has a given constant modulus in advance compatible with the characteristics of the motor on cubesat B. So we pose: c being a given constant such that: m(t).c<10 −4 N. We suppose that the mass of the cubesat B is equal to 3 kg. Then c N k g ,  < --As we will use, in our numerical application, the kilometre as distance unit, we can take c 0.3 10 km s 0.3 10 ms .
Hence, the constant c has the acceleration dimension. It is the thrust acceleration. The equations of motion are: for any u(t)). We suppose that B is on an elliptic orbit. At the instant t 1 , we suppose that B is on the xy plane; but this initial position will be modified after the adjustment of orbits, in order to achieve the rendezvous of B with A.
The choice of the direction of the thrust vector T  will be specified according to the desired result. Let us define a local reference frame in B, r h , , , q (ˆˆˆ) such that:

Adjustment of angular momentum vectors
To obtain the coincidence of the two angular momentum vectors, in an optimal way, we pose h h and we calculate the derivative of the square of the modulus of this vector with respect to time, where the scalar product of two vectors is denoted by a dot: Consequently, the optimal thrust to make h B

  coincide with h A
  is to be collinear to vector U  and in the opposite direction. (a) Let us calculate the derivative of the eccentricity vector for a given thrust.
We have, omitting the index B and writing μ for μ T : inally, we get: Consequently, for a thrust T  of given modulus, the maximum decrease of e D    | | is obtained when the thrust is collinear to the vector U '

Adjustment of the two orbits
It is possible to choose an optimal direction of thrust to make e B  and h B   coincide respectively with e A  and h .

 
For that, we will consider the sum: We introduce the positive factor r (expressed in km) such that Let us calculate the derivative of S r with respect to time:   are given by (4) and (6) So, we have: negative with maximum modulus, at each point of the orbit. We expect that, in the beginning of the convergence process, h D    | |will decrease up to be near , e D    | | then S r will tends rapidly to zero. We make the calculation in Fortran (FT95) to test the convergence speed, for the adjustment of both orbits, for different values of .
r We notice that the number t j of days necessary to obtain this adjustment is a function of r which has a unique minimum for some r = r 1 . Hence, the optimal direction of the thrust is obtained for ρ=ρ 1 .
Finally, for ρ=ρ 1 the direction of T  becomes optimal. Then, the integration of differential equations (E), insures, in a minimal time, the convergence of S e h 2 2 | | }to zero and the coincidence of B orbit with A orbit.

Numerical results
We recall that the mass of the Earth is equal to M T =5.9736×10 24 kg, the constant of the gravitation is equal to G= 6.67×10 −11 N·m 2 ·kg −2 , the geocentric gravitational constant is equal to μ T =G·M T =398 600.4418 (km) 3 s −2 . The radius of the Earth is equal to 6378 km.
We suppose that B is, at the instant t 1 , on the xy plane; but this position will be modified after the adjustment of orbits, in order to achieve the rendezvous of B with A.
We take the example where B has a circular orbit: e B =0, r B =6982 km, i B =69°and Ω B =237°. We suppose that B is, at the instant t 1 (calculated below), in the xy plane where r B  makes with Ox axis an angle equal to Ω B =237°. The coordinates of B, for tt 1 , can be written before the instant t 1 of the thrust beginning: x t r cos cos n t t cos i sin sin n t t y t r sin cos n t t cos i cos sin n t t z t r sin i sin n . t t where n r .
( ) The method that we propose, in this paper, is valid for any choice of the satellite A. We will carry out the numerical applications by taking as satellite A, the satellite ALMASAT-1 whose orbit parameters are: where i A is the inclination of the orbit of A, Ω A is the angle made by the line of nodes with the x-axis of the Galilean frame centred on the Earth, ω A is the argument of the latitude of the pericentre counted from the line of the ascending nodes and M A is the mean anomaly which is written: where t p is the instant passage of ALMASAT-1 to pericentre and a A is the semi-major axis of the orbit of A. We know that the ALMASAT- We take, in the following, t p =0, so we have: t 1 =3277.931 913 s≈3277.932 s. Moreover, having μ T and n A , we obtain the semi-major axis a A : T a n 4 7252.394 km.
then we obtain, at the instant t 1 , the value α A1 of the angle α=ω A +θ: w 2.962403336 0.7862045057 2.176198830 rad, We deduce also the value r A1 of r A : r a 1 e 1 e . cos a 1 e cos E 7816.537634 km.

Calculation program
We specify that the calculation program performs the following operations: α) It integrates the differential equations by the Runge-Kutta method of the fourth order, every second, from time t 1 , over the interval [t 1 +i−1, t 1 +i], i=1, 2, ... . q where ω B is the argument of the latitude of B pericentre. α A =ω A +q A , where ω A is the argument of the latitude of A pericentre.

Adjustment of the two orbits
We saw, above, that the optimal direction of the thrust is such that: u c , For a given , r we will make the numerical calculation from the instant t 1 , with this choice of T  which leads to the adjustment of the two orbits. So, we take, at the initial instant t 1 : [  β) The calculation is stopped when the following conditions are verified: 10 km s and 10 10 Then we notice that the time t j , necessary to obtain the adjustment of the two orbits, has, as a function of , r a unique minimum near ρ 2 =1000. For this value of ρ, the numerical calculation shows that the orbit of B coincide with the orbit of A at the instant t 2 =18 144 389 s = 210 days 0 h 6 min 29 s.

Achievement of the rendezvous α)
We recall that, for ρ=ρ 1 =1150, we make the calculation using a FORTRAN 95 program, with initial conditions at the instant t 1 , M A1 =191°.986 and M B1 =0. We stop the calculation when the two conditions (9) are verified. Under these conditions, the B orbit coincide with the A orbit, for t 2 =17 982 828 s ≡ 208 days 3 h 13 m 48 s, namely 208 days 2 h 19 m 10 s after the instant t 1 . But at the instant t 2 , we obtain M B2 =184°.816, whereas the mean anomaly of A is: M A2 =239°.690.
Thus a final step is necessary to achieve the rendezvous: we must choose a convenient position of B on his initial circular orbit, to make B coincide with A at the final time, t 2 .
The satellite B is behind the satellite A when B arrives at the instant t 2 in the orbit of A, with the following difference between mean anomalies: Consequently, at the instant t , 1 the mean anomaly of B must be: M B1 =58°.092 798. We expect, with this new initial condition that the cubesat B will coincide with A at the instant t 2 . But the calculation, made with this new initial condition for the mean anomaly of B, shows a drift for t 2   We stop the calculation when the conditions (9) are verified. By repeating this calculation, we verify that |δ M B | decreases rapidly and tends to zero. Finally, when we take M B1 =52°.181 928, the rendezvous is obtained, at t 2 =17 981 421 s ≡ 208 days 2 h 50 m 21 s, where M A2 -M B2 =0°.000 05, with a distance between A and B approximately equal to: 7252×(0.000 05×π/180) km≈6.32 m.
We summarize this result in the following table 1, where the convergence of this process is obtained after 9 steps: β) Let us use the same algorithm when the calculation is stopped under the conditions (10). We take ρ =ρ 2 =1000, and the initial conditions: M A1 =191°.986 and M B1 =0, at the instant t 1 . Then the orbit of B coincide with the orbit of A, for t 2 =18 144 389 s ≡ 210 days 0 h 6 min 29 s, namely 209 days 23 h 11 min 51 s after the instant t 1 . We obtain at the instant t 2 : We repeat this calculation until δM B1 tends to zero. The convergence of this process is obtained after 8 steps, as we can see in the following Remark: To achieve the rendezvous, after the adjustment of orbits, an alternative to this procedure is to use the method of local approach exposed in our previous work 8 .

Conclusion
We present in this paper a new method which permits to a cubesat, equipped with a low-thrust propulsion system which covers an orbit around the Earth, to make rendezvous with a given satellite on an elliptic orbit around the Earth. The thrust is very low but continuous; its acceleration is constant and its direction is adjusted continuously to obtain optimal convergence of the angular momentum vector and the eccentricity vector of the cubesat respectively to the angular vector and the eccentricity vector of the satellite. We gave a numerical application of this method taking ALMASAT-1 as a satellite, where the cubesat starts on circular orbit. Evidently this method can be applied for any satellite around the Earth or another planet. It can be used to remove old satellites or debris around the Earth by deorbiting them, in the aim to burn the debris in the dense atmosphere.

Appendix A. Some formulas
The eccentricity vector can be written, on the local frame as follows, where θ is the true anomaly: We have in the Galilean frame, where the index B is omitted to simplify the notations: We take ρ=1150 and we stop the calculation when conditions (9) are verified 9 . Then the orbit of B coincides with the orbit of A at the instant t 2 =17 981 421 s, namely 208 days 2 h 50 m 43 s after the instant t 1 of the beginning of the thrust. The following curves show that a B , e B , i B , Ω B and ω B tend respectively to a A , e A , i A , Ω A and ω A ; they coincide for t=t 2 .
1) Semi-major axis a B :      We take ρ=1000. The calculation is stopped when the conditions (10) are verified 10 . The rendezvous is obtained at t 2 =18 142 340 s, or 209 days 22 h 37 m 42 s after the instant t 1 . The following curves show that a B , e B , i B , Ω B , ω B and M A tend respectively to a A , e A , i A , Ω A , ω A and M A ; they coincide for t=t 2 .
1) The semi-major axis a B of B: 2) The eccentricity of B: 3) The inclination i B of B: Figure C1. The achievement of the rendezvous. The semi-major axis a B as a function of t, for ρ=1000. Figure C2. The achievement of the rendezvous. The eccentricity e B as a function of t, for ρ=1000. Figure C3. The achievement of the rendezvous. The inclination i B as a function of t, for ρ=1000.