Session 6: Analytical Approximations for Low Thrust Maneuvers

Session 6 : Analytical Approximations for Low Thrust Maneuvers As mentioned in the previous lecture, solving non-Keplerian problems in general requires the use of perturbation methods and many are only solvable through numerical integration. However, there are a few examples of low-thrust space propulsion maneuvers for which we can find approximate analytical expressions. In this lecture, we explore a couple of these maneuvers, both of which are useful because of their precision and practical value: i) climb or descent from a circular orbit with continuous thrust, and ii) in-orbit repositioning, or walking. Spiral Climb/Descent We start by writing the equations of motion in polar coordinates, d2 r −r dt2 � dθ dt �2 + µ = ar r2 (1) d2 θ 2 dr dθ aθ + = 2 r dt dt r dt (2) We assume continuous thrust in the angular direction, therefore ar = 0. If the acceleration force along θ is small, then we can safely assume the orbit will remain nearly circular and the semi-major axis will be just slightly different after one orbital period. Of course, small and slightly are vague words. To make the analysis rigorous, we need to be more precise. Let us say that for this approximation to be valid, the angular acceleration has to be much smaller than the corresponding centrifugal or gravitational forces (the last two terms in the LHS of Eq. (1)) and that the radial acceleration (the first term in the LHS in the same equation) is negligible. Given these assumptions, from Eq. (1), dθ ≈ dt � µ r3 d2 θ 3 ≈− 2 dt 2 → � µ dr r5 dt (3) Substituting into Eq. (2), we obtain a differential equation for r, which can be integrated directly for an initial radius r0 and time t0 , r r0 dr ≈ r3/2 t t0 2aθ √ dt µ 1 → 1/2 r0 − 1 r1/2 aθ ≈ √ (t − t0 ) µ (4) Re-arranging the integrated expression and setting t0 = 0, r ≈ r0 (1 − aθ t/v0 )2 where the velocity of the initial circular orbit is v0 = also notice that, t µ/r0 . From the definition of Δv, we � aθ dt = aθ t ≈ Δv = 0 1 (5) µ − r0 � µ r (6) From Eq. (5) we observe that the trajectory will be a climbing, or descending spiral, de­ pending on whether aθ is positive or negative. Eq. (6) shows that the Δv defined in terms of the perturbation acceleration is equal to the change in velocity between the initial and final orbits. To note that this change of velocity is not equal to the ideal rocket Δv, we compare it with the corresponding Hohmann transfer, 2µr − r0 (r + r0 ) ΔvH = (vp − v0 ) + (v − va ) = � r µ r0 r � + µ − r 2µr0 r(r + r0 ) (7) By definition, this is an impulsive maneuver, not susceptible to losses from gravitational, or any other externally applied force. In addition, these impulses are provided at both ends of the apsidal line resulting in the optimal Δv to change the altitude and circularize the orbit. Assume the final and initial orbital radii are related by r = nr0 and calculate the ratio Δv/ΔvH , Δv = ΔvH � " �−1  �  √ � 2 n 2 1+ −1 n+1 (8) This ratio is, as expected, always larger than unity, meaning that Hohmann transfers always require a lower Δv. Eq. (8) is shown in the figure below. The logarithmic plot is symmetric about n = 1, accounting for descending (n < 1) and ascending (n > 1) trajectories. 2.2 2 ∆V/∆VH 1.8 1.6 1.4 1.2 1 −1 10 0 10 1 10 n = r/r0 2 10 The spiral approximation holds as long as the orbit remains near-circular, and within our assumptions this will be true if the angular acceleration is small. We need to be careful when applying these results, since in many instances we are interested in constant acceler­ ation maneuvers for which the strength of the angular acceleration relative to gravity and centrifugal forces will decrease for n < 1, and increase for n > 1. Therefore, descending spirals could safely be analyzed with these tools whereas the approximation to an ascending spiral will eventually fail. 2 In any event, low thrust spiral maneuvers are not optimal in the sense that work is always done against the gravitational field. In this particular case, the thrust vector is not perfectly aligned along the “circular” trajectory and a small eccentricity will be introduced, which will increase with time. It is interesting to note that the same result in Eq. (5) could be obtained with alternative methods. For instance, considering again that the orbit remains near circular, we can calculate the rate of change of the orbital energy, and make this equal to the thrust power delivered to the vehicle, dET F~ · ~v = dt m d  µ − = aθ dt 2r → r µ r (9) Eq. (9) is identical to Eq. (4) and therefore its solution is the same. The spiral trajectory appears to be a trivial solution, but there are some subtleties. Notice that the velocity increment ∆v is actually equal to the decrease in orbital velocity. The rocket is pushing forward, but the velocity is decreasing. This is because in a r−2 force field, the kinetic energy is equal in magnitude but opposite in sign to the total energy, 1 2 µ µ µ µ v = ET + = − + = 2 r 2r r 2r (10) Eq. (5) suggests that, in principle, escape conditions will be reached at t = v0 /aθ when r → ∞. But of course, the orbit is no longer near-circular when approaching escape, so we cannot expect this result to be precise. We could obtain a more precise determination of escape conditions (∆vesc ) in the following way. The radial velocity ṙ can be calculated from Eq. (5) by differentiation. Notice that this is in the nature of an iteration, since r̈ was neglected in the energy balance which led to Eq. (9). We then obtain, ṙ = 2aθ r0 /v0 (1 − aθ t/v0 )3 (11) The angular component vθ = rθ˙ is approximately the orbital velocity, i.e., rθ˙ = r µ = r r   µ aθ t − aθ t = v0 1 − r0 v0 (12) The overall kinetic energy is therefore,  v2 1 2 v2 = ṙ + r2 θ˙2 = 0 2 2 2 " aθ t 1− v0 2 2 4 (aθ r0 /v02 ) + (1 − aθ t/v0 )6 # (13) Escape conditions are reached when the total energy vanishes, i.e., 1 2 µ v − = 0 2 r 1 2 or 3  v v0 2 − r0 = 0 r (14) Substituting, 1 2  aθ tesc 1− v0 2  2 2 2 (aθ r0 /v02 ) aθ tesc + − 1− = 0 v0 (1 − aθ tesc /v0 )6 so we have, aθ tesc 1− = v0  2aθ r0 v02 1/4  = 2aθ µ/r02 1/4 = (2ε)1/4 (15) where ε = aθ /(µ/r02 ) is the ratio of thrust to gravitational accelerations, and as before should be small for the approximation to hold. Since ∆vesc = aθ tesc ,   ∆vesc ≈ v0 1 − (2ε)1/4 (16) This result is useful in obtaining a preliminary determination of escape conditions, but because of the assumptions and the eventual increase of the eccentricity, we do not expect that Eq. (16) will converge to the exact result, even for very small values of ε. To evaluate this model, the equations of motion are solved numerically and the total energy is tracked until it vanishes. At that point, we compute the quantities shown in the table below. ε dr ds
esc ∆vesc v0 sesc r0 resc √ ε r0 1−∆vesc /v0 ε1/4 10−2 0.5327 0.7615 51.13 0.8518 0.7541 10−3 0.5346 0.8657 503.6 0.8535 0.7552 10−4 0.5348 0.9245 5011.4 0.8538 0.7549 10−5 0.5347 0.9575 50036 0.8534 0.7554 In the table, s is the distance along the orbital trajectory. In particular, ~s · F~ is the work done to reach escape and should be equal to the orbit’s initial energy, sesc = v2 µ/2r0 = 0 aθ 2aθ (17) The factor 21/4 = 1.19 in Eq. (16) is definitely larger than the values in the last column of the table above. In consequence, a better expression to use would be,   ∆vesc ≈ v0 1 − 0.754ε1/4 (18) Finally, we observe that the escape radius will be given by, 0.85 resc ≈ r0 √ ε and the rate of climb compared to the distance traveled at escape is, 4 (19) � dr ds � ≈ 0.53 (20) esc which is very far from a circular orbit, and closer to a trajectory that starts to move away from the orbit focus in almost a straight line. The numerical results would be slightly different if instead of using angular thrust we use tangential thrust. The analysis, however, is not as straightforward and closed analytical expressions like those shown above are more difficult to obtain. Re-positioning in Orbits: Walking Suppose now that we want to move a satellite in a circular orbit to an angular position Δθ apart in the same orbit, in a time Δt (assumed to be several orbital times at least). The general approach is to transfer to a lower (for Δθ > 0) or higher (for Δθ < 0) nearby orbit, then drift in this faster (or slower) orbit for a certain time, then return to the original orbit. The analysis is similar for low and high thrust, because in either case the satellite is nearly in the same orbit even during thrusting periods, and as we found out for spiral transfers, the Δv for orbit transfer is equal to the magnitude of the difference between the beginning and ending orbital speeds. In detail, of course, if done at high thrust the maneuver involves a two-impulse Hohmann transfer to the drift orbit and one other two-impulse Hohmann transfer back to the original orbit. For the low-thrust case, continuous thrusting is used during both legs, with some guidance required to remove the very slight radial component of iv picked up during spiral flight (and ignored here). We will do the analysis for the low-thrust case only, then adapt the result for high-thrust. Let δθ be the advance angle relative to a hypothetical satellite remaining in the original orbit and left undisturbed. The general shape of the maneuver is sketched below: forward thrust ing st coa backwards thrust The orbital angular velocity is Ω = µ/r3 , and its variation with orbit radius is, 3 dr d(δθ) dΩ = − Ω = 2 r dt The radial variation can be computed through the power balance, d � µ � Fi · iv − = = ac dt 2r m � µ r 5 → 1 dr = 2ac r dt (21) � r µ (22) The rate of change of the angular velocity is,  r  dΩ d2 (δθ) 3 r0 3ac = = − Ω 2a = − c dt dt2 2 µ r0 (23) in which we have made r ≈ r0 as an approximation, since we do not expect the radius to change significantly during the maneuver, d(δθ) 3ac = − t+A dt r0 (24) where A is a constant. Starting from t = 0, δθ = 0, then d(δθ)/dt = 0, we have, d(δθ) 3ac =− t dt r0 → δθ = − 3ac 2 t 2r0 (t < t1 ) (25)   t1 t− 2 (26) After t = t1 , we continue to drift at a constant rate, d(δθ) 3ac =− t1 dt r0 and since we start from, δθ(t1 ) = − 3ac 2 t 2r0 1 the relative angle δθ during the coasting phase is, δθcoast 3ac 2 3ac 3ac t1 = − t1 − t1 (t − t1 ) = − r0 2r0 r0 At the end of coasting , and we have, 3ac t1 δ θ(∆t − t1 ) = − r0  3t1 ∆t − 2  (27) and, after a second period t1 of reversed thrust, we return to the initial orbit with d(δθ)/dt = 0, and with δθ as in Eq. (27), plus a further δθ(t1 ). The total ∆θ is, 3ac t1 ∆θ = − r0   3t1 3ac t21 ∆t − − 2 2r0 → ∆θ = − 3ac t1 (∆t − t1 ) r0 (28) Clearly, the mission (given ∆θ and ∆t) can be accomplished with different choices of thrusting time t1 (but notice that t1 < ∆t/2 in any case). The required acceleration ac and ∆v = 2|ac |t1 depend on this choice, ac = − r0 ∆θ 3t1 (∆t − t1 ) 6 (29) Δv = 2 r0 Δθ 3 (Δt − t1 ) (30) Not surprisingly, we find again that low thrust ends up as a penalty on Δv, so that the thrusting time should be selected as short as possible within the available on-board power. In the limit of impulsive thrust, we realize that Δt cannot really be any less than the Hohmann transfer time. A more detailed analysis of this case confirms that, for the high thrust case, Eqs. (29-30) are indeed valid with t1 = π/Ω. The power per unit mass required is, P 1 |F |c r0 Δθ c = = m 2η m 6ηt1 (Δt − t1 ) (31) Finally, some analyses might benefit from expressing the results in terms of the coasting time tc = Δt − 2t1 , so that, t1 = Δt − tc , 2 Δt − t1 = Δt + tc 2 and t1 (Δt − t1 ) = Δt2 − t2c 4 We then have, ac = − 4r0 Δθ 3(Δt2 − t2c ) (32) Δv = 4 r0 Δθ 3 (Δt + tc ) (33) P 2r0 Δθ c = m 3η(Δt2 + t2c ) (34) coasting reduces Δv, but increases P/m (not much if tc /Δt is small). 7