16.851 Problem Set #2 (Solutions)

TO: Prof. John Keesee, Prof. David Miller, and Ms. Marilyn Good SUBJECT: Problem Set #2 (Orbit and Propulsion) Problem Set #2 (Orbit and Propulsion Subsystems) Motivation In a mission design in which an orbit transfer is necessary, the orbit transfer and propulsion system types may drive of the design. Many missions desire minimal mass to reduce the mission cost, which may be achieved in part by choosing the appropriate propulsion system. Some military missions and manned missions, however, may prefer minimal orbit transfer time over minimal cost associated with minimal mass. In many missions, the choices for the orbit transfer and propulsion system may be obvious. For some missions, the appropriate design is not obvious, given that designing a propulsion system depends on the transfer orbit chosen, and vice versa. For such design problem, a tool that designs the orbit transfer and propulsion system may alleviate the difficulty in solving the coupled problem. Problem Statement The objective is to design a tool that compares a set of feasible orbit transfer and propulsion system combinations for a given mission requirement. The mission requirement is in terms of the desired orbit transfers. For simplicity, we will only consider an orbit transfer from a lower circular orbit to a higher circular orbit. The output of the design tool should be a set of feasible transfer orbit and propulsion system combinations. The comparison metric is the mass of the propulsion system and the time required for orbit transfer. Note that the set of possible transfer orbits depends on the propulsion system type, and the total ∆Vtotal required depends on the transfer orbit type. While this problem can become very complex depending on the accuracy of the computation and analysis, a right set of assumptions should be applied to minimize the computational complexity while maintaining the fidelity for comparison among various designs (e.g. assume fixed dry mass, but excluding the propulsion system mass). The satellite mission will have the following guidelines: 1. Orbit transfer takes place purely by provision of additional energy by the propulsion system. No plane change is considered. 2. Shapes of initial and final orbits are circular. 3. Both initial and final orbits are assumed to satisfy all components which affect possibility, duration, and efficacy of satellite mission, such as launch windows, radiation effects, earth convergence and others. 4. The dry mass, mdry, (excluding the propulsion system mass) is fixed. The power subsystem mass have a fixed mass, although its crucial influence on comparison between electric propulsion systems are noted. 5. The maximum initial mass, mmax, is fixed. 6. Transfer orbit types to consider are: a. Elliptical & Hyperbolic Transfer b. Hohmann Transfer c. Hohmann Transfer Segments d. Spiral Transfer 7. Types of propulsion system to consider are: a. Solid Motor b. Chemical i. Bipropellant ii. Monopropellant c. Electric Propulsion systems can be broken down further with different propellants and associated Isps. cold gas systems are excluded because of the combination of low thrust and low Isp. In general, cold gas systems are not applicable to orbit transfers even for satellites that are extreme sensitivity to contamination and/or for which the complexity of the propulsion system is of an issue. Approach Flow chart below shows the outline of out approach procedure. The correlation between orbit and propulsion system is complementary for each other. The velocity change calculated from initial and final orbits selects propulsion systems such that achieve that value. On the other hand, propulsion system with a certain velocity change budget constrains the types of orbital maneuvers. Calculation starting with certain assumptions and releasing those assumptions as variable later is required to get out of this loop relationship. In this MATLAB code, we commenced the calculation with propulsion system and associated achievable velocity change limit as inputs, and checked its validity with velocity change obtained from transfer orbits later. Propulsion System Type Mass Assumptions: Maximum total mass, Fixed dry mass Propulsion Calculate Maximum impulsive ∆V achievable ∆V max Initial Orbit Final Orbit ∆V max Decide applicable Transfer Orbit Types Given Orbit Hohmann Transfer Elliptical & Hyperbolic Transfer Hohmann Transfer Segments Spiral Transfer Transfer Time, and total ∆V Calculate required propulsion mass to achieve ∆V Iterate Propulsion Propulsion mass including Tank STEP 1: Compute the maximum ∆V available. To calculate the maximum velocity change achievable and the maximum firing time by a propulsion system, Impulsive Maneuvers for Orbit Transfer is considered. For most orbit transfer calculations, each ∆V is assumed impulsive given the thrust is “relatively” high. That is, with high thrust capability, an infinitesimal time is necessary to achieve the necessary ∆V, and as such, we can approximate the burn as an impulse. If the thrust is relatively low, however, the maneuver will require a finite amount of time. Figure 1 illustrates the case in which the burn time is finite. In such case, the thrust vector is not along the flight path and results in gravity loss as shown by the second term on the right hand side of Eq. (1). If the thrust vector were aligned along the flight path, the gravity loss would be eliminated. The resulting orbit transfer, however, will not be the intended one. Thus, for a spacecraft to achieve the intended orbit transfer, the burn time must be short enough such that the gravity loss is negligible compared to the thrust. ∆V = gI sp ln tf mi − ∫ g sin γ dt ti mf (1) θ mg ti burn time flight path tf γ T Figure 1. Forces during finite burn. We use a small angle approximation, and assume that for γ ≤ π/12, i.e. θ ≤ π/6, the gravity loss effect is minimal. Then, given the gravitational constant µ of the central body, current circular orbital radius r, mass of the spacecraft mi, and available propellant mass mprop, we can compute the maximum achievable ∆Vmax: ⎛ ⎞ mi ∆Vmax = gI sp ln ⎜ ⎟ ⎝ mi − ∆mmax ⎠ (2) ⎡ θ Tr 3/ 2 ⎤ ∆mmax = min ⎢ m p , ⎥. gI sp µ ⎥⎦ ⎢⎣ (3) where Note that mo and mp are initially unknown. Thus, we assume mi = mmax and mp = mmax – mdry. These values must be update once the actual values are known, and the remaining design processes must be repeated. This is implemented in compute_impulsive_Delta_V_max.m. Input r T Isp initial total mass, i.e. mdry + mwet available propellant mass gravitational constant of the body orbiting radius of the initial orbit thrust specific impulse ∆Vmax maximum impulsive ∆V available mi mprop µ Output STEP2: Compute the ∆V required for the orbit transfer. Calculate the velocity change required and transfer time for each orbit transfer from orbit elements of given initial and final orbit elements. Verify that velocity change for each firing is smaller than the maximum velocity achievable by a propulsion system. Orbit element values provided by given initial and final orbits are follows: • rA: Radius of initial orbit • rB; Radius of final orbit • ∆Vmax; Maximum velocity change according to propellant The satellite velocities on initial, transfer, and final orbits are calculated using these values mentioned above and the gravitational constant of a body the spacecraft is orbiting. Hohmann Transfer Hohmann Transfer orbit’s ellipse is tangent to both the initial and final orbits at the transfer orbit’s perigee and apogee respectively. This is one of the most fuel-efficient transfers between two circular coplanar orbits and the orbit transfer time is relatively small. The disadvantage is that a failure is unrecoverable during the orbit transfer. Atx: Semimajor axis of transfer ellipse atx = rA + rB 2 Via: Velocity in the initial orbit viA = Vfb: Velocity in the final orbit v fB = µ rA µ rB ⎛2 1 ⎞ Vtxa: Velocity at perigee point on transfer orbit vtxA = µ ⎜ − ⎟ ⎝ rA atx ⎠ Vfxb: Velocity at apogee point on transfer orbit ⎛2 1 ⎞ vtxB = µ ⎜ − ⎟ ⎝ rB atx ⎠ The velocity changes required at each staging is computed as follows. ∆VA: Velocity change at the perigee of transfer orbit (1st burn) ∆VA = VtxA − ViA ∆VB: Velocity change at the apogee of transfer orbit (2nd burn) After verifying that both ∆VA and ∆VB are smaller than transfer time can be determined by: ∆V max from ∆VB = V fB − VtxB propellant, total ∆V and ∆Vtotal = ∆VA + ∆VB ∆t = π atx 3 µ This is implemented in compute_Hohmann_transfer.m. Input rA rB µ inner orbit radius outer orbit radius gravitational constant of the body orbiting ∆V ∆t ∆V required for transfer time required for transfer Output Elliptical & Hyperbolic Transfer Elliptical & Hyperbolic Transfer is a modified version of Hohmann to achieve rapid transfer. Extra energy was put at the first burn, and transfer orbit can be a bigger ellipse than Hohmann transfer orbit or a hyperbola. First, decide maximum rc, apoapsis of transfer orbit, which should achieve the fastest transfer, according to ∆V max due to propellant. atx; Semimajor axis of transfer eclipse atx = Via: Velocity in the initial orbit Via = Vfb: Velocity in the final orbit Vfb = ra + rc 2 µ ra µ rb 1 ⎞ ⎛ 2 Vtxa: Velocity at perigee point on transfer orbit Vtxa = µ ⎜ − ⎟ ⎝ ra atx ⎠ Vfxb: Velocity at apogee point on transfer orbit 1 ⎞ ⎛2 Vtxb = µ ⎜ − ⎟ ⎝ rb atx ⎠ ∆Va : Velocity change at the perigee of transfer orbit (1st burn) ∆Va = Vtxa − Via ∆Vb : Velocity change at the apogee of transfer orbit (2nd burn) ∆Vb = Vfb 2 + Vtxb 2 − 2VfbVtxb cos φ where ⎡ e sin ν ⎤ ⎣1 + e cosν ⎥⎦ φ = tan −1 ⎢ ( ) ⎡ ⎧ atx 1 − e 2 ⎫⎤ − 1⎬ ⎥ ⎢⎨ rb ⎭⎥ ν = cos −1 ⎢ ⎩ ⎢ ⎥ e ⎢ ⎥ ⎢⎣ ⎥⎦ Considering that both ∆Va and ∆Vb are smaller than ∆V max from propellant, we can get maximum rc, accordingly. Also, e is also determined by: e = 1− ra atx 1) ∆Va > ∆Vb First, the code assumes that ∆Va is the limitation. By deploying the equation of ∆Va = ∆V max , maximum atx and rc are calculated as follows: 1 atx = 2 ⎛ µ ⎞ + ∆Vmax ⎟ ⎜ 2 ⎝ rA ⎠ − µ rA rc = 2atx − ra Verify that ∆Vb corresponding to this rc is smaller than ∆V max . 2) ∆Vb > ∆Va Try out this case only if ∆Vb > ∆V max . Start from the decent number of rc (>rb), and increase the number until ∆Vb value becomes bigger than ∆V max . The maximum value of rc which still achieve ∆Vb < ∆V max determines the fastest transfer. In either case, plugging into rc value obtained into the equations above leads to ∆V total and transfer time required. ∆Vtotal = ∆Va + ∆Vb 3 t = 0.001583913atx 2 (E − e sin E ) where ⎡ e + cosν ⎤ E = cos −1 ⎢ ⎣1 + e cosν ⎥⎦ This is implemented in compute_high_energy_transfer.m. Input µ ∆Vmax inner orbit radius outer orbit radius gravitational constant of the body orbiting maximum impulsive ∆V available ∆V ∆t ∆V required for transfer time required for transfer rA rB Output Hohmann Transfer Segments This transfer orbit is often chosen for a low-thrust transfer, using low thrust chemical or electrical propulsion. First, calculate the value of atk (semimajor axis of k th transfer eclipse) achievable using ∆V max from propellant by deploying ∆Vak = ∆V max : atk: Semimajor axis of k th transfer 1 atk = 2 ⎛ ⎛ 2 ⎞ ⎞ 1 ⎜ µ⎜ − ⎟ + ∆V max ⎟ ⎜ ⎜⎝ ra at (k − 1) ⎟⎠ ⎟ 2 ⎝ ⎠ − µ ra rck: Apoapsis of k th transfer orbit rck = 2atx − ra where at0 = ra (k=0) atK = rb (k=K, final orbit) Via: Velocity in the initial orbit Via = Vfb: Velocity in the final orbit Vfb = µ ra µ rb 1 ⎞ ⎛ 2 Vtka: Velocity at perigee point on k th transfer orbit Vtka = µ ⎜ − ⎟ ⎝ ra atkx ⎠ ∆Vi Velocity change at the perigee of transfer orbit on k th burn ∆Vak = Vtka − Vt (k − 1) a Total velocity change and total transfer time required are accumulated in each transfer orbit segments, and are determined by: k 1 ⎞ ⎛ 2 ∆Vtotal = ∑ ∆Vi + µ ⎜ − ⎟ ⎝ rb atk ⎠ i =1 k −1 t = ∑ 2Π (ra + rk )3 8µ i =1 (ra + rk ) 3 +Π 8µ Spiral Orbit Transfer Given a propulsion system with a very low acceleration capability (i.e.T/W < 10-3), we cannot approximate each maneuver as impulsive due to high gravity loss (see Approach, STEP1). This implies that the aforementioned orbit transfers cannot be used. In such case, we can use a spiral orbit transfer (see Figure 2). During a spiral orbit transfer, the spacecraft thrusts continuously until it reaches the desired orbit. rf ro Figure 2. Low thrust spiral orbit. Edelbaum [2] has introduced a method for optimal low-thrust transfer between inclined circular orbits, which is also discussed in [1]. For coplanar orbit transfer between two circular orbits, however, we can derive the equations for the necessary ∆V and the time for orbit transfer using basic knowledge of the orbital mechanics and few key assumptions. Assuming that the spacecraft thrusts in the direction of it’s velocity vector, the change in the specific mechanical energy, ε, of the spacecraft with mass m is dε T = ⋅v = a⋅v , dt m (4) where T is thrust (assumed constant), a is the acceleration of the spacecraft, and v is velocity. Given that the thrust is very small, we can assume that the orbit remains relatively circular. That is, dε µ dr = 2 , dt 2r dt (5) where µ is the gravitational constant and r is the radius at which the spacecraft is in. Equating Eqs. (4) and (5), and assuming constant acceleration, we can compute the amount of time the spacecraft takes, ∆t, to transfer to a circular orbit at rB radius: ∆t = 1⎛ µ µ⎞ − ⎜⎜ ⎟ a ⎝ rA rB ⎟⎠ (6) The total ∆V necessary to transfer to a circular orbit with radius rB from rA is then given by: ∆Vtotal = a ⋅ ∆t (7) This is implemented in compute_spriral_transfer.m. Input mi T inner orbit radius outer orbit radius gravitational constant of the body orbiting initial total mass, i.e. mdry + mwet thrust ∆V ∆t ∆V required for transfer time required for transfer rA rB µ Output STEP 3: Compute the propulsion system mass. Once propulsion system type has been chosen and the necessary ∆V budget is determined, we can compute the necessary propellant mass mprop from the rocket equation: m prop = mdry e ∆V g0 I sp − mdry (8) where Isp is the specific impulse of the propulsion system chosen, g0 is the acceleration due to Earth’s gravitational force at the sea level. Note that in general, the true dry mass mdry is unknown since the propulsion system mass is unknown among others. Thus, this number has to be estimated in the beginning and must be updated as better estimate of the dry mass is generated. Now with the propellant mass, we can calculate the mass of the propulsion system. Once the propulsion system mass is calculated, the dry mass must be updated and the whole process must be iterated until it converges to a desired level of confidence. Chemical We could go through a detailed analysis to determine the mass of a chemical propulsion system, such as monopropellant and bipropellant systems. The analysis, however, become very complex and tedious. Furthermore, there are many detailed design choices required to compute the propulsion system mass. While such tool would be great for future use, we take a very simple approach to approximating the chemical propulsion system mass. That is, we assume that approximately 85 ~ 90% of the system mass is the propellant mass [5]. Pressurant, tanks, lines, fittings, components, etc. contributes to the remaining 15%. mTotal = m prop (9) 0.85 Solid Similar to the chemical propulsion system, we take a crude approximation of the solid rocket propulsion system mass as a mass fraction. Typically, 82 ~ 94% of a solid rocket prolusion system mass account for the propellant mass. We take the conservative 82%. mTotal = m prop (10) 0.82 Electric Propulsion For electric propulsion, we will use a simple empirical model [3] to estimate the mass of the propulsion system. This is will more than suffice for our desired level of fidelity. This tool could be replaced with higher fidelity model as desired or necessary in the future. First, we compute the efficiency of the propulsion system η as follows: η = ηthrusterη PowerProcessor = A + B ln ( I sp ) (11) where the constants A and B are defined in Table 1. Table 1. Empirical data for modeling the efficiency and the specific mass of an electric propulsion system. Propulsion System H2 Arcjet NH3 Arcjet Ar Ion A -2.024 Constants for Models B C 5.0 1.8 0.307 4490 D 0 0 -0.781 Xe Ion Hg Ion Ar MPD Ar PIT -1.776 -0.765 -0.591 -1.99 0.307 0.181 0.126 0.32 123100 82870 7 7 -1.198 -1.136 0 0 Given the flow rate of the propellant as a function of the propellant mass and the duration for with the thruster is turned on ∆tburn, m = m prop (12) ∆tburn we can compute the necessary power: P= m ( I sp g 0 ) 2 (13) 2η Then, the mass of the thruster and the power processor is computed from an empirical formula using Table 1 where: mThruster + mPowerProcessor = CI sp D P (14) Finally the tank mass is calculated using one of the equation listed in Table 2. Table 2. Electric propulsion system propellant tank mass. Propellant Type NH3 NH3 He Xe Propellant Mass Range (kg) 5000 ~ 18300 18300 ~ 2200 5000 ~ 13000 5000 ~ 2200 Tank Mass (kg) 120 + 0.173 mprop + 2.28 mprop 2/3 1020 + 0.198 mprop 610 + 0.493 mprop 52 + 0.075 mprop + 0.154 mprop2/3 The total electrical propulsion system ass is: mTotal = mThruster + mPr ower Pr ocessr + mTank + m prop (15) While the mass of the necessary power source for the propulsion system may be significant, it isn’t included as it overlaps with the power subsystem’s requirements. Implementation Source Code compute_impulsive_Delta_V_max.m function Delta_V_max = compute_impulsive_Delta_V_max(m_o,m_prop,mu,r,thruster) % DELTA_V_MAX = COMPUTE_IMPULSIVE_DELTA_V_MAX(M_O,M_PROP) % % Inputs: % M_O Initial total mass, i.e. m_dry + m_wet % M_PROP Available propellant mass % MU Gravitational constant of the body orbiting % R radius of the initial orbit % TRUSTER Thruster info % % Outputs: % DELTA_V_MAX Maximum impulsive Delta V available % reject any invalid inputs if (m_o <= m_prop) error('Error! Propellant mass must be less than the initial spacecraft mass.'); end theta = 2*pi/12; % Angle of the orbit for wich the Delta V occurs. [rad] % Assume pi/12 ~ 0, i.e. small angle approximation. g = 98.1; % Earth gravity at sea level [m/s^2] % Maximum propellant available during an impulsive burn. [kg] Delta_m_max = min(m_prop, theta * thruster.T * r^(3/2) /... (g*thruster.I_sp * sqrt(mu))); % Maximum Delta_V available during an impulsive burn. [m/s] Delta_V_max = g * thruster.I_sp * log(m_o/(m_o - Delta_m_max)); compute_high_energy_transfer.m function [Delta_V,Delta_t] = compute_high_energy_transfer(mu,r_A,r_B,Delta_V_max) % [DELTA_V,DELTA_T] = COMPUTE_HIGH_ENERGY_TRANSFER(MU,R_A,R_B,Delta_V_max) % % Inputs: % MU Gravitational constant of the body orbiting % R_A Radius of the inner circular orbit % R_B Radius of the outer circular orbit % Delta_V_max Maximum impulsive Delta_V available % % Outputs: % DELTA_V Total Delta_V necessary for Hohmann transfer % Delta_T Orbit transfer time % Constrain so that the first burn is no greater than the max allowed. Delta_V_A = Delta_V_max; v_iA = sqrt(mu/r_A); v_fB = sqrt(mu/r_B); % Velocity for inner circular orbit [m/s] % Velocity for outer circular orbit [m/s] % Compute the transfer orbit assuming that the maximum Delta_V_max is used % at point A. v_txA = Delta_V_max + v_iA; % Velocity for transfer orbit at first burn [m/s] a_tx = 1/(2/r_A - v_txA^2/mu); % Semi-major axis of the transfer orbit [m] e = 1 - r_A/a_tx; % Eccentricity of the transfer orbit nu = acos((a_tx*(1-e^2)/r_B - 1)/e);% True anomaly at the second burn [rad] phi = atan(e*sin(nu)/(1+e*cos(nu)));% Flight path angle at second burn [rad] v_txB = sqrt(mu*(2/r_B - 1/a_tx)); % Velocity for transfer orbit at second burn [m/s] % Second burn Delta_V necessary [m/s] Delta_V_B = sqrt(v_fB^2 + v_txB^2 - 2*v_fB*v_txB*cos(phi)); % Check to see if the second burn is less than the maximum allowed. if (Delta_V_B <= Delta_V_max) % The first burn is the limiting factor. Delta_V = Delta_V_A + Delta_V_B; % Total Delta_V [m/s] E = acos((e + cos(nu))/(1+e*cos(nu))); % Eccentricy anomaly at second burn [rad] Delta_t = sqrt(a_tx^3/mu)*(E-e*sin(E)); % Orbit transfer time [sec] return else % The second burn is the liminting factor. Delta_V_B = Delta_V_max; % Delta_V at second burn [m/s] a_tx = search_a_tx(mu,r_A,r_B,Delta_V_B); % Semi-major axis of transfer oribt [m] v_txA = sqrt(mu*(2/r_A - 1/a_tx)); % Velocity for tranfer orbit at first burn [m/s] Delta_V_A = v_txA -v_iA; % Delta_V at first burn [m/s] Delta_V = Delta_V_A + Delta_V_B; % Total Delta_V [m/s] e = 1 - r_A/a_tx; % Eccentricity of the transfer orbit nu = acos((a_tx*(1-e^2)/r_B - 1)/e); % True anomaly at the second burn [rad] E = acos((e + cos(nu))/(1+e*cos(nu))); % Eccentricy anomaly at second burn [rad] Delta_t = sqrt(a_tx^3/mu)*(E-e*sin(E)); % Orbit transfer time [sec] end function Delta_V_B = compute_Delta_V_B(a_tx,mu,r_A,r_B) e = 1 - r_A/a_tx % Eccentricity of the transfer orbit nu = acos((a_tx*(1-e^2)/r_B - 1)/e); % True anomaly at the second burn [rad] phi = atan(e*sin(nu)/(1+e*cos(nu))); % Flight path angle at second burn [rad] v_txB = sqrt(mu*(2/r_B - 1/a_tx)); % Velocity for transfer orbit at second burn [m/s] v_fB = sqrt(mu/r_B); % Velocity for outer circular orbit [m/s] Delta_V_B = sqrt(v_fB^2 + v_txB^2 - 2*v_fB*v_txB*cos(phi)) function a_tx = search_a_tx(mu,r_A,r_B,Delta_V_B) allowable_percent_error = 0.000001; a_tx_low = r_B/2; a_tx = (r_A + r_B)/2; a_tx_high = (r_A + r_B)*10000; quit = 0; while (not(quit)) difference = Delta_V_B - compute_Delta_V_B(a_tx,mu,r_A,r_B); % This should be zero if ((a_tx_high - a_tx)/a_tx < allowable_percent_error | ... (a_tx - a_tx_low)/a_tx < allowable_percenT_error); quit = 1; elseif (difference > 0) a_tx_low = a_tx; a_tx = (a_tx_high - a_tx)/2; elseif (difference < 0) a_tx_high = a_tx; a_tx = (a_tx - a_tx_low)/2; end end compute_Hohmann_transfer function [Delta_V,Delta_t] = compute_Hohmann_transfer(mu,r_A,r_B) % [DELTA_V,DELTA_T] = COMPUTE_HOHMANN_TRANSFER(MU,R_A,R_B) % % Inputs: % MU Gravitational constant of the body orbiting % R_A Radius of the inner circular orbit % R_B Radius of the outer circular orbit % % Outputs: % DELTA_V Total Delta_V necessary for Hohmann transfer % Delta_T Hohman transfer time % reject any invalid inputs if (r_A > r_B) error('Error! Inner orbit radius is larger than the outer.'); end v_iA = sqrt(mu/r_A); v_fB = sqrt(mu/r_B); a_tx = (r_A + r_B)/2; v_txA = sqrt(mu*(2/r_A - 1/a_tx)); v_txB = sqrt(mu*(2/r_B - 1/a_tx)); Delta_V_A = v_txA - v_iA; Delta_V_B = v_fB - v_txB; Delta_V = Delta_V_A + Delta_V_B; Delta_t = pi*sqrt(a_tx^3/mu); % % % % % % % % % Velocity for inner orbit [m/s] Velocity for outer orbit [m/s] Semi-major axis of transfer orbit [m/s] Velocity at periapsis for m/s] Velocity at apoapsis for [m/s] Delta_V at periapsis of [m/s] Delta_V at apoapsis of transfer [m/s] Total Delta V required transfer [m/s] Hohmann transfer time [s] compute_spiral_transfer function [Delta_V,Delta_t] = compute_spiral_transfer(mu,r_A,r_B,m_i,T) % [DELTA_V,DELTA_T] = COMPUTE_SPIRAL_TRANSFER(MU,R_A,R_B,M_I,T) % % Inputs: % MU Gravitational constant of the body orbiting % R_A Radius of the inner circular orbit % R_B Radius of the outer circular orbit % M_I Initial Total Mass % T Thrust % % Outputs: % DELTA_V % Delta_T Total Delta_V necessary for Hohmann transfer Hohman transfer time % reject any invalid inputs if (r_A > r_B) error('Error! Inner orbit radius is larger than the outer.'); end if (T/m_i*g > 10e3) error('Error! Too much thrust for spiral transfer assumptions.'); end a = T/m_i; % Acceleration [m/s^2] Delta_t = (sqrt(mu/r_A) - sqrt(mu/r_B))/a; % Orbit transfer time [s] Delta_V = a*Delta_t; % Total Delta V required for orbit transfer [m/s] References 1. Chobotov V., Orbital Mechanics, 2nd ed. Reston, VA: American Institute of Aeronautics and Astronautics, Inc., 1996. 2. Edelbaum, T., “Propulsion Requirements for Controllable Satellites,” ARS Journal, Vol. 31, August 1961, pp. 1079-1089. 3. Humble R., Henry G., Larson W., Space Propulsion Analysis and Design., New York: The McGraw-Hill Companies, Inc. Primis custom Publishing, 1995. 4. Sutton G. and Biblarz O., Rocket Propulsion Elements, New York: John Wiley & Sons, Inc., 2001. 5. Wertz, J., and Larson, W., Space Mission Analysis and Design, 3rd ed. Torrance, CA: Microcosm Press, 1999.