Next Article in Journal
On the Wake Properties of Segmented Trailing Edge Extensions
Next Article in Special Issue
Leading-Edge Roughness Affecting Diamond-Wing Aerodynamic Characteristics
Previous Article in Journal
Cost-Effectiveness of Structural Health Monitoring in Fuselage Maintenance of the Civil Aviation Industry
Previous Article in Special Issue
Flight Load Assessment for Light Aircraft Landing Trajectories in Windy Atmosphere and Near Wind Farms
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

On Four New Methods of Analytical Calculation of Rocket Trajectories

by
Luís M. B. C. Campos
*,† and
Paulo J. S. Gil
CCTAE, IDMEC, Instituto Superior Técnico, Universidade de Lisboa, Av. Rovisco Pais, 1049-001 Lisboa, Portugal
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Aerospace 2018, 5(3), 88; https://doi.org/10.3390/aerospace5030088
Submission received: 15 June 2018 / Revised: 20 July 2018 / Accepted: 27 July 2018 / Published: 15 August 2018
(This article belongs to the Special Issue Feature Papers in Aerospace)

Abstract

:
The calculation of rocket trajectories is most often performed using purely numerical methods that account for all relevant parameters and provide the required results. There is a complementary need for analytical methods that make more explicit the effect of the various rocket and atmospheric parameters of the trajectory and can be used as test cases with unlimited accuracy. The available analytical methods take into account (i) variable rocket mass due to propellant consumption. The present paper includes four new analytical methods taking into account besides (i) also (ii) nonlinear aerodynamic forces proportional to the square of the velocity and (iii) exponential dependence of the mass density with altitude for an isothermal atmospheric layer. The four new methods can be used in “hybrid analytical-numerical” approach in which: (i) the atmosphere is divided into isothermal rather than homogeneous layers for greater physical fidelity; and (ii) in each layer, an exact analytical solution of the equations of motion with greater mathematical accuracy than a numerical approximation is used. This should allow a more accurate calculation of rocket trajectories while discretizing the atmosphere into a smaller number of layers. The paper therefore concentrates on four analytical methods of calculation of rocket trajectories in an isothermal atmospheric layers using new exact solutions of the equations of motion beyond those currently available in the literature. The four methods are developed first for the simpler case of a vertical climb and will be subsequently extended to the practically more relevant case of a gravity turn.

1. Introduction

Four new methods of analytical calculation of rocket trajectories are presented. The analytical methods have at least three advantages. First, the equations show explicitly the influence on the rocket trajectory of all parameters in the model. Second, the analytical solutions can be calculated with unlimited accuracy, within the domain of validity of the model. Thirdly, in the case of the present models, the numerical computation of analytical results can be performed in a few seconds on an ordinary PC, and thus computational cost is almost negligible. The numerical methods are irreplaceable for more complex models with higher computational cost and the interpretation of results is not always clear. Thus, analytical and numerical methods are complementary, with the former providing “exact” test cases for, and helping the interpretation, of the latter. Analytical and numerical methods have common features; for example, they can be extended from a single stage to a multi-stage rocket by taking the final conditions of the preceding stage as the initial conditions for the next stage. The numerical methods that dominate the research literature on rocket trajectories account for up to thirteen effects listed below in the next paragraph. The analytical methods that appear mostly in textbooks include three or four effects. The models in the present paper account for seven effects as indicated in Table 1 and thus help to bridge the gap between the textbooks and research literature.
When considering the trajectory of rockets in the atmosphere, there are at least thirteen effects which may be of interest: (i) the variation of rocket mass with time, due to propellant consumption, which depends on the throttle schedule, and is most important if the propellant fraction of initial mass is large (it can be up to 90%); (ii) the angle of the thrust with the rocket axis, which may be zero, small or non-negligible, in particular in the case of thrust vectoring, when it can vary with time; (iii) the variation of thrust over time, not only due to throttling, but also due to variation of gas expansion characteristics in the exhaust nozzle, as external pressure decreases with altitude; (iv) the aerodynamic drag, that can significantly reduce acceleration and speed in the atmosphere at lower altitudes; (v) the aerodynamic lift, which is most significant for winged transatmospheric vehicles, and can generate significant forces normal to the flight path; (vi) the angle-of-attack, which by tilting the aerodynamic axes relative to the flight path, causes the lift and drag to have both longitudinal and normal force components; (vii) the dependence of aerodynamic forces, viz. lift and drag (and also side force), on velocity, which is nonlinear, viz. quadratic, in the simplest subsonic approximation; (viii) the dependence of aerodynamic forces on atmospheric mass density which varies significantly with altitude; (ix) the effect of wind that varies with geographical location and time of day; (x) non-uniform gravity, if the altitude range of flight is not small compared to the Earth radius; (xi) the Coriolis force which can cause deviation from plane trajectories for large velocities; (xii) the effect of the centrifugal force which is significant for large distance and angular velocity relative to the center of the Earth; and (xiii) the discarded mass at each stage of a multi-stage rocket, leading to a succession of flight phases, each taking as initial conditions the final conditions of the preceding, with the mass reduced by that of the discarded stage.
Of the preceding thirteen effects, for atmospheric flight over an altitude range much smaller than the radius of the Earth h R = 6378   k m , it is usual to omit (x, xi, xii), that is take uniform gravity and neglect centrifugal and Coriolis forces. The effect (xiii) of multiple rocket stages reduces to a sequence of one-stage problems, which poses new challenges mainly for trajectory optimization. The remaining effects (i) to (ix) are taken into account in the present calculation of rocket trajectories, requiring for definiteness that some assumptions be made: (i) the propellant flow rate c is assumed to be constant, so that the rocket mass decreases linearly with time from the initial m 0 to the residual m 1 mass; (ii) the thrust is assumed to make a constant angle with the rocket axis, usually small; (iii) the variation of thrust with external pressure is neglected, so it is taken as constant; (iv–vi) lift (iv) and drag (v) are considered, assuming (vi) a constant angle of attack relative to the flight path; (vii) the aerodynamic forces, viz. lift and drag, are assumed to be proportional to the square of the rocket velocity and to the atmospheric mass density; (viii) the mass density is assumed to decay exponentially with altitude, corresponding to an isothermal atmospheric layer. The effect of wind speed (ix) can be included in the aerodynamic forces. The ensemble of ten effects included in the analytic models of rocket trajectories presented in this paper is compared in Table 1 with the three or four effects considered in the literature [1,2,3,4,5], mostly consisting of textbooks.
The set of assumptions made in the present paper for the analytic calculation of rocket trajectories is far less restrictive than those available in the preceding literature. For example, the standard textbook solutions concern only the effect of mass variation due to propellant consumption, and do not take into account the dependence of the aerodynamic forces on airspeed and altitude. These and other effects can be readily included in numerical methods for engineering purposes. The analytical methods are complementary in that they can provide more insight into the effects of the various rocket, aerodynamic, and atmospheric parameters on the trajectory. Thus, this paper reduces the large gap between the numerical methods and the available analytical solutions whose domain of application is currently much more restricted. The more general analytical methods can also be used in hybrid analytical-numerical approaches, using analytical solutions in a smaller number of atmospheric layers than a purely numerical approach. The four novel analytical methods to calculate rocket trajectories use series expansions in terms of (I) rocket mass, (II) propellant mass, (III) time, or (IV) altitude, with details provided next.
The trajectories of rockets in the atmosphere are specified by a set of ordinary differential equations with time as the independent variable, coupling the path coordinates (or angles) and their first and second-order derivatives, viz. velocities and accelerations, translational for coordinates and angular for angles. The coupled system of ordinary differential equations for the trajectories is: (i) nonlinear if the aerodynamic forces are included because they are quadratic functions of the velocity, in the simplest case; (ii) have coefficients depending on time if propellant consumption is taken into account; and (iii) have coefficients depending on position if the atmospheric mass density in the aerodynamic forces is taken as a function of altitude. Thus, under simple but realistic conditions, the trajectories of rockets are specified by coupled nonlinear ordinary differential equations, with variable coefficients depending on position and time; this renders analytical solutions difficult and requires some care in numerical algorithms, in order to avoid singularities and ensure convergence. A possibly effective approach is a semi-analytical method, using simple analytical solutions as inputs to more complex numerical algorithms.
In the absence of aerodynamic or propulsive side force, and omitting the Coriolis force, the rocket trajectory lies on a plane. Two choices of coordinate systems are the most common: (a) a Cartesian earth-fixed coordinate frame; (b) an orthogonal Frenet–Serret frame moving along the trajectory with a longitudinal coordinate tangent to it, and hence the transversal coordinate towards or away from the center of curvature. Using the moving frame (b) and taking as variables the velocity and flight path angle, the acceleration components along and across the flight path are nonlinear; however, in this case, it is straightforward to eliminate time from the trajectory by dividing the longitudinal by the transversal acceleration, leading to a single differential equation for the path. This method was originally applied [6] to calculate the ballistic trajectory of a projectile of constant mass in a uniform gravity field, with drag proportional to a power of the velocity and opposite to it; the extension of this solution to include the effects of thrust and variable mass has been researched extensively [7,8]. The simplicity of this method is lost if the dependence of aerodynamic forces on altitude through the mass density is included; in this case, altitude appears explicitly in the equations of the trajectory, requiring an additional integral equation, without which its solution is not closed. The original method (a) can still be applied by altitude “steps” with constant density.
The present paper aims at including the continuous dependence of aerodynamic forces on altitude through the mass density, and chooses as a starting point the dynamical trajectory equations (Section 2), using an earth-fixed Cartesian frame (Section 2.3) rather than the orthogonal frame moving along the trajectory (Section 2.1); in the latter case, altitude can be eliminated (Section 2.2) by introducing the radius of the curvature of the trajectory. The use of fixed rather than a moving reference frame leads to two rather than one differential equation for the trajectory (Section 3). The problem is reduced to a single differential equation (Section 3.1) by the assumption of constant flight path angle. This rather restrictive assumption is used to develop more simply a one-dimensional method of solution (Section 3.2), specifying the altitude as a function of time (Section 3.3). A particular case of trajectory with a constant flight path angle is the vertical climb, which is quite inefficient at gaining speed because it works entirely against gravity; for this reason, the gravity turn is used [1,2,3,4,5] which could be approximated by a sequence of straight paths, each with a constant flight path angle. This restriction will be removed in subsequent work, where a continuous gravity turn is considered, using an extended variant of the methods developed here, in order to determine the curved two-dimensional trajectory, by specifying the horizontal and vertical coordinates and velocities as functions of time.
Two methods of solution of the equations of motion are considered, both using the atmospheric mass density as altitude variable; the time variable is: (Section 3) the burned propellant mass fraction in method I; and (Section 4) the residual mass fraction in method II. Both methods involve the same dimensionless thrust, gravity and drag parameters (Section 3.1 and Section 4.1), and use series expansions (Section 3.2 and Section 4.2), whose leading terms are calculated explicitly (Section 3.3 and Section 4.3). It is shown that method II is (Section 5.1) more accurate for a long time, e.g., burn-out, and method I for a short time, e.g., initial climb. Since method I has a more tedious algebra than method II, an alternative method III most suited for a short time is presented (Section 5.2) using a Taylor series of altitude as a function of time. All three methods apply to all stages of a multi-stage rocket, since they allow arbitrary initial altitude and velocity. The case of free flight without thrust post burnout is considered by a fourth method IV, extending the known solution for constant atmospheric mass density [3], to exponential variation of altitude; this method leads to a series expansion with the terms of all orders calculated explicitly (Section 5.3).

2. Dynamical Equations of the Rocket Trajectory

The dynamical equations of the rocket trajectory are written in two reference frames: (i) an earth-fixed Cartesian frame (Section 2.3) for which the components of acceleration and gravity are simplest; (ii) a moving Frenet–Serret orthogonal frame with one axis tangent to the trajectory at rocket velocity. In the latter case (ii), two possibilities are considered, viz: (ii-a) eliminating time (Section 2.1), in which case the altitude remains in the mass density in the aerodynamic forces; (ii-b) eliminating the atmospheric mass density (Section 2.2) by using the lift-to-drag ratio, in which case the radius of curvature of the trajectory appears.

2.1. Orthogonal Frame Moving with the Velocity

The geometry of a rocket trajectory (Figure 1) may be used to establish the equations of motion in an orthogonal frame, with one axis tangent to the trajectory, moving at the rocket velocity. The tangential Equation (1a) and normal Equation (1b) force balance are specified [9,10] respectively by:
m d V d t = m g sin γ D cos α L sin α + T cos ( α + ε ) ,
m V d γ d t = m g cos γ D sin α + L cos α + T sin ( α + ε ) ,
where: (i) on the right-hand side, the inertia force is the product of the mass by the acceleration, whose tangential component involves only the total velocity V, whereas the normal component involves also the flight path angle γ ; (ii) on the right-hand side. The first term is the gravity force, which is vertical downwards, and thus makes an angle π / 2 γ with the velocity V, so that cos ( π / 2 γ ) = sin γ appears in Equation (1a) and sin ( π / 2 γ ) = cos γ appears in Equation (1b); (ii) the second term is the drag force D, which makes an angle π + α with the velocity, where α is the angle-of-attack, so that cos ( π + α ) = cos α appears in Equation (1a) and sin ( π + α ) = sin α appears in Equation (1b); (iii) the third term is the lift force L, which is orthogonal to the drag, and makes an angle π / 2 + α with the velocity, so that cos ( π / 2 + α ) = sin α appears in Equation (1a) and sin ( π / 2 + α ) = cos α appears in Equation (1b); (iv) the last term is the thrust T, which makes an angle ε with the rocket axis, hence an angle α + ε with the velocity, so that cos ( α + ε ) appears in Equation (1a) and sin ( α + ε ) in Equation (1b).
Of the coefficients in the equations of motion: (i) the mass m varies with time due to propellant consumption, e.g., for a constant flow rate is given by
c d m d t = const . : m ( t ) = m 0 c ( t t 0 ) ,
where the initial m 0 and final m 1 mass determine the burn time t 1 ,
t 0 t t 1 : t 1 t 0 = m 0 m 1 c ,
after which the trajectory is unpowered; (ii) this implies that the thrust T is zero for t > t 1 , and before it is constant and non-zero if the effect of atmospheric pressure on nozzle expansion ratio is neglected, otherwise it is a function of altitude T ( z ) ; (iii) the acceleration of gravity may be taken as constant for altitudes much smaller than the Earth radius z 1 z 0 R = 6378   k m ; and (iv) the remaining forces in Equations (1a,b) are lift and drag:
L = 1 2 ρ ( z ) S ( V u ) 2 C L ( α ) ,
D = 1 2 ρ ( z ) S ( V u ) 2 C D ( α ) ,
where the dimensionless lift C L and drag C D coefficients depend on the angle of attack α but not the cross-section S; in the simplest case, if weakly affected by the Mach and Reynolds numbers, they are independent of velocity V and mass density ρ , which thus appear respectively quadratically and linearly in the aerodynamic forces Equations (3a,b); the wind speed u projected along the tangent to the flight path is subtracted from the velocity V in an Earth fixed reference frame to specify the airspeed V u in the aerodynamic forces, namely lift Equation (3a) and drag Equation (3b); the atmospheric mass density decays rapidly with altitude, viz. exponentially with scale height
ρ 0 ρ ( z 0 ) : ρ ( z ) = ρ 0 exp z z 0 ,
for an isothermal atmospheric layer.
Substituting Equation (4) in Equations (3a,b), the equations of motion (1a,b) become, with the preceding assumptions:
m 0 c ( t t 0 ) d V d t = ( m 0 c t ) g sin γ + T ( z ) cos ( α + ε ) 1 2 S ( V u ) 2 ρ 0 e ( z z 0 ) / C D ( α ) cos α + C L ( α ) sin α ,
m 0 c ( t t 0 ) V d γ d t = ( m 0 c t ) g cos γ + T ( z ) cos ( α + ε ) + 1 2 S ( V u ) 2 ρ 0 e ( z z 0 ) / C L ( α ) cos α C D ( α ) sin α ,
which are: (i) nonlinear because of the normal acceleration and aerodynamic forces; (ii) have coefficients depending on time through the mass varying with propellant consumption, and also possibly thrust vectoring ε ( t ) and angle-of-attack variation α ( t ) for flight control; (iii) the coefficients depend on altitude, through the atmospheric mass density in the aerodynamic forces, and nozzle pressure ratio in the thrust; and (iv) the wind speed is a local weather variable depending on time or altitude along the flight path.
If the angle-of-attack is zero α = 0 , the thrust is aligned with the rocket axis ε = 0 , and lift is neglected compared to drag L D , then the equations of motion (1a,b) simplify to:
α = 0 = ε , L D : m d V d t = m g sin γ + T D ,
m V d γ d t = m g cos γ ,
which have been considered extensively in the literature [1,2,3,4,5]; eliminating d t by taking the ratio leads to
1 V d V d γ = tan γ + D m g sec γ T m g sec γ ,
which is a nonlinear differential equation relating the velocity V and flight path angle γ . The latter can be generalized by taking the ratio of Equations (1a,b), viz.
1 V d V d γ = T cos ( α + ε ) m g sin γ D cos α L sin α T sin ( α + ε ) m g cos γ D sin α + L cos α .
Even if the angle-of-attack α ( t ) and thrust vector angle ε ( t ) are constant in Equation (8), or zero in Equation (7), the time dependence remains in the mass m ( t ) , e.g., Equations (2a,b) for constant rate of propellant consumption. The altitude dependence remains in Equation (8) through the thrust being a function of nozzle pressure ratio and aerodynamic forces depending on mass density; even for constant thrust and zero lift in Equation (7), the altitude still appears through the atmospheric mass density Equation (4) in the drag force Equation (3a). The altitude z is related to the velocity and flight path angle via an integration over time,
z z 0 = t 0 t V ( t ) sin γ ( t ) d t .
The inclusion of Equation (9) in the mass density Equation (4) through the aerodynamic forces Equations (3a,b), turn Equation (7) or Equation (8) into an integro-differential equation for V , γ . This could be avoided by approximating Equation (7) or Equation (8) as a sequence of altitude steps with constant density, along the rocket trajectory (Figure 2).

2.2. Influence of the Lift-to-Drag Ratio

The altitude dependence could be eliminated from the coefficients of the general equations of motion (1a,b) in the form:
D cos α + L sin α = T cos ( α + ε ) m g sin γ m d V d t ,
D sin α L cos α = T sin ( α + ε ) m g cos γ m d γ d t ,
by dividing to obtain
λ + tan α λ tan α 1 = T cos ( α + ε ) m g sin γ m d V d t T sin ( α + ε ) m g cos γ m d γ d t ,
where the aerodynamic forces appear only through the drag-to-lift ratio
λ D L = C D ( α ) C L ( α ) = C D 0 C L ( α ) + K ¯ + K C L ( α ) ,
the latter is independent of velocity and wind speed, and also of mass density and hence also altitude. Thus, for constant thrust, the altitude is eliminated from Equation (11), and the drag-to-lift ratio Equation (12) depends only on angle-of-attack, e.g., in the simplest case it involves, besides the lift coefficient C L ( α ) the: (i) friction drag coefficient C D 0 ; (ii) the lift-induced drag coefficient K; and (iii) the coefficient K ¯ for quadratic non-parabolic lift-drag polar. Writing the drag-to-lift ratio in terms of an aerodynamic angle δ in Equation (13a):
λ D L = tan δ ,
λ + tan α λ   tan α 1 = tan ( α + δ ) ,
and the left-hand side of Equation (11) takes the form Equation (13b). Figure 1 shows the rocket axis making an angle-of-attack α with the velocity vector V , and the lift L and drag D components of the total aerodynamic force F ; the aerodynamic angle δ is the angle of the total aerodynamic force with the lift, and appears to be added to the angle-of-attack in Equation (13b), the left-hand side of the equation of the trajectory (11).
On the right-hand side of Equation (11), the differential of time d t can be eliminated by using the arc length d s and the radius of curvature Equation (14a) of the trajectory viz.:
R = d s d γ ,
V d γ d t = V d γ d s d s d t = V 2 R ,
in the normal Equation (14b) and tangential Equation (14c) acceleration,
d V d t = d V d γ d γ d s d s d t = V R d V d γ .
Substituting Equations (14b,c) in Equation (11) leads to an equation for the trajectory,
tan ( α + δ ) = λ + tan α λ tan α 1 = T cos ( α + ε ) m g sin γ m V R d V d γ T sin ( α + ε ) m g cos γ m V 2 R ,
which also involves the velocity V and flight path angle γ like Equation (8) but eliminates dependence on altitude, and quadratic terms on the velocity and wind speed in the aerodynamic forces.
Using the simplifying assumptions as in Equations (6a,b) of zero angle-of-attack and thrust vector angle, and negligible lift, simplifies Equation (15) to
α = 0 = ε , L D : tan ( α + δ ) = λ + tan α λ tan α 1 = V d V d γ + g R sin γ R T m V 2 + g R cos γ ,
which should be compared with Equation (7), just as the more general forms, respectively Equations (15) and (8), can be compared: (i) both are first-order nonlinear differential equations for the velocity V and flight path angle γ ; (ii) the coefficients depend on time in both because of variable mass m ( t ) in Equations (7) and (16), and possibly also thrust vectoring ε ( t ) and angle-of-attack scheduling α ( t ) in Equations (8) and (15).
The differences are that: (i) the altitude appears explicitly in Equations (7) and (8), even at constant thrust, leading by Equation (9) to an integro-differential equation; (ii) the altitude does not appear in Equations (15) or (16) for constant thrust, but the radius of curvature Equation (14a) appears, viz.
R = d s d γ = d s d t d γ d t = V γ ˙ ,
the radius of curvature is the ratio of the velocity to the angular velocity of the flight path angle. The radius of curvature is a function of the arc length R ( s ) , which is related to the velocity in the absence of wind by
s = 0 t V ( t ) 1 d t ,
and substitution of Equation (18) in Equation (15) or Equation (16) leads an to integro-differential equation. This can be avoided by approximating the trajectory by circular arcs, which have constant radius of curvature (Figure 3); the flight path angle can be made continuous at the junction of two circular arcs, by choosing the center at the intersection of the normals at the end-points.

2.3. Fixed Cartesian Frame with Axis Along Altitude

In an earth-fixed Cartesian frame with a z-axis along the altitude and the x-axis down range, the acceleration is linear in the equations of motion:
m x ¨ = D cos ( γ + α ) L sin ( γ + α ) + T cos ( γ + α + ε ) ,
m z ¨ = m g D sin ( γ + α ) + L cos ( γ + α ) + T sin ( γ + α + ε ) ,
and the acceleration of gravity appears as a constant and only in one equation; the remaining forces, viz. thrust, lift and drag are similar in Equations (1a,b) and (19a,b), with the flight path angle γ added in all trigonometric coefficients. The aerodynamic forces Equations (3a,b) can be written Equations (20a,b) in the absence of wind [11,12]:
L = 1 2 ρ ( z ) S ( x ˙ 2 + z ˙ 2 ) C L ( α ) ,
D = 1 2 ρ ( z ) S ( x ˙ 2 + z ˙ 2 ) C D ( α ) ,
and using the same assumptions Equations (2a,b) and (4) as before, the equations of motion (19a,b) become:
m 0 c ( t t 0 ) x ¨ = T cos ( γ + α + ε ) 1 2 ρ 0 S x ˙ 2 + z ˙ 2 e ( z z 0 ) / × C D ( α ) cos ( γ + α ) + C L ( α ) sin ( γ + α ) ,
m 0 c ( t t 0 ) ( z ¨ + g ) = T sin ( γ + α + ε ) 1 2 ρ 0 S x ˙ 2 + z ˙ 2 e ( z z 0 ) / × C D ( α ) sin ( γ + α ) C L ( α ) cos ( γ + α ) ,
where the flight path angle
tan γ = d z d x = z ˙ x ˙
must be substituted.
The substitution of the flight path angle in the equations of motion (21a,b) is made using:
cos γ = x ˙ 2 + z ˙ 2 1 / 2 x ˙ ,
sin γ = x ˙ 2 + z ˙ 2 1 / 2 z ˙ ,
in:
x ˙ 2 + z ˙ 2 1 / 2 cos ( γ + α ) = x ˙ cos α z ˙ sin α ,
x ˙ 2 + z ˙ 2 1 / 2 sin ( γ + α ) = x ˙ sin α + z ˙ cos α ,
and likewise for cos , sin ( γ + α + ε ) replacing α by α + ε . Thus, the equations of motion (21a,b) in Cartesian coordinates Equations (24a,b) become
[ m 0 c ( t t 0 ) ] x ˙ 2 + z ˙ 2 1 / 2 x ¨ z ¨ + g = T ( z ) cos ( α + ε ) sin ( α + ε ) sin ( α + ε ) cos ( α + ε ) x ˙ z ˙ 1 2 ρ 0 S exp [ ( z z 0 ) / ] x ˙ 2 + z ˙ 2 C D cos α + C L sin α C D sin α + C L cos α C D sin α C L cos α C D cos α + C L sin α x ˙ z ˙ .
The time dependence appears in the mass, angle-of-attack α and thrust vector angle ε ; the altitude dependence appears in the atmospheric mass density and thrust.
The trajectory can be approximated by at least three alternative methods: (i) altitude steps (Figure 2) of constant atmospheric mass density (Section 2.1); (ii) circular arcs (Figure 3) of constant radius of curvature (Section 2.2); and (iii) straight segments with a constant flight path angle (Figure 4). This approximation will be released in subsequent work, and is used next to develop more simply a method of trajectory calculation. The assumption of constant flight path angle γ in Equation (22) implies that the trajectory is specified by one equation only, e.g., Equation (21b) for altitude:
m 0 c ( t t 0 ) ( z ¨ + g ) = T ( z ) sin ( γ + α + ε ) + 1 2 ρ 0 S z ˙ 2 csc 2 γ exp [ ( z z 0 ) / ] × C L ( α ) cos ( γ + α ) C D ( α ) sin ( γ + α ) ,
which allows for vertical flight γ = π / 2 but excludes horizontal flight γ = 0 . The opposite would be the case choosing Equation (21a).
The flight path angle can be taken as constant only if the vertical and horizontal accelerations are also in the ratio Equation (22), viz.
x ¨ sin γ = z ¨ cos γ .
Substituting Equations (19a,b) in Equation (27) leads to the identity
m g cos γ + D sin α = L cos α + T sin ( α + ε ) ,
which states (Figure 1) the balance of forces normal to the flight path. In particular, vertical flight γ = π / 2 or cos γ = 0 is possible for zero angle of attack α = 0 , thrust aligned with the rocket axis ε = 0 and negligible lift L = 0 . If the angle of attack is zero α = 0 and lift is neglected L = 0 , then thrust at angle ε 0 to the rocket axis would lead to a trajectory at an angle π / 2 γ to the vertical, with flight path angle
α = 0 = L : cos γ = T m g sin ε ,
which would be constant for constant mass and fixed thrust angle, or mass variation compensated by thrust vectoring. In general, the aerodynamic forces vary with velocity and altitude in Equation (28), and rocket mass changes with time, so that a constant flight path angle could be maintained by thrust modulation or vectoring or angle-of-attack scheduling. In this case, Equation (28) could be used together with the equation of motion tangent to the trajectory,
m z ¨ csc γ = m g sin γ + T sin ( α + ε ) + L cos α D sin α ,
or alternatively:
m 0 c ( t t 0 ) z ¨ csc γ = m g sin γ + T ( z ) sin ( α + ε ) + 1 2 ρ 0 S z ˙ 2 csc 2 γ exp [ ( z z 0 ) / ] C L ( α ) cos α C D ( α ) sin α ,
with the assumptions Equations (2a), (3a,b), and (4).

3. Trajectory Calculation via a Power Series of Time

The method of trajectory calculation using power series of the time is presented most simply for a vertical climb, which involves only three dimensionless parameters representing gravity, thrust and drag (Section 3.1); the ratio of atmospheric mass densities at arbitrary altitude z relative to initial altitude z 0 is used as a dimensionless altitude variable, and the time is made dimensionless using the initial mass and propellant rate (Section 3.2); the nonlinear differential equation and initial conditions determine the leading coefficients of a power series solution.

3.1. Dimensionless Gravity, Thrust and Drag Parameters

For zero angle-of-attack α = 0 , neglecting lift L = 0 , and thrust aligned with the rocket axis, the condition Equation (28) of a straight flight path is met for cos γ = 0 or γ = π / 2 , i.e., a vertical climb, in which case Equation (30a) simplifies to
α = ε = 0 = L : m z ¨ = m g + T D ,
which balances the inertia force against gravity, thrust and drag, leading in the absence of wind to
m 0 c ( t t 0 ) ( z ¨ + g ) = T 1 2 C D S ρ 0 z ˙ 2 exp [ ( z z 0 ) / ] ,
using the assumptions Equations (2a,3b,4). The initial altitude is taken to be z 0 at time t = t 0 so that the results can be applied to any stage of a multistage rocket. Thus, the differential equation is solved with the initial conditions:
z ( t 0 ) = z 0 ,
z ˙ ( t 0 ) = v 0 ,
specifying altitude Equation (33a) and velocity Equation (33b) at time t = t 0 of thrust initiation for an upper stage, or lift-off for the first stage.
Next, dimensionless independent and dependent variables are introduced. The dependent variable is altitude, z, which is taken relative to initial altitude, made dimensionless Equation (34a) by dividing by the atmospheric scale height:
ζ ( z z 0 ) / ,
τ m ( t ) m 0 1 = c ( t t 0 ) m 0 ,
whereas time is replaced as the independent variable by Equation (34b) the ratio of the mass of propellant burned c ( t t 0 ) to the initial mass m 0 . Substituting Equations (34a,b) in the equation of motion (32), it takes the dimensionless form
( 1 τ ) ζ + a = b f ζ 2 e ζ ,
where prime denotes derivative with regard to τ , viz. ζ = d ζ d τ , and the three parameters:
a m 0 2 g c 2 ,
b m 0 T c 2 ,
f C D ρ 0 S 2 m 0 ,
are also all dimensionless.
The dimensionless parameters specify the relative importance of gravity Equation (36a), thrust Equation (36b), and drag Equation (36c), which are in the proportion
m o g : T : C D 2 ρ 0 S m 0 c 2 m 0 .
The last term in Equation (37) involves the reference mass m * defined by Equation (38a) as the mass of a column of air, with the mass density ρ 0 at the initial altitude z 0 , height equal to the atmospheric scale height , and cross-section S equal to that of the rocket:
m * = ρ 0 S ,
f = C D 2 m * m 0 ,
the dimensionless drag parameter Equation (36c) equals Equation (38b) the reference mass normalized to the initial mass m * / m 0 , multiplied by half the drag coefficient. The dimensionless acceleration Equation (36a) is the ratio Equation (38c) of the acceleration of gravity g to the reference acceleration Equation (38d):
a g a * ,
a * = c 2 m 0 2 ,
which is an acceleration constructed from the propellant rate Equation (2a), initial m 0 and atmospheric scale height in Equation (4). The dimensionless thrust parameter Equation (36b) is similar to Equation (36a) replacing the initial rocket weight m 0 g bt the thrust T.

3.2. Atmospheric Mass Density as a Function of the Mass of Burned Propellant

In order to eliminate the exponential nonlinearity from the differential equation of motion (35), the initial atmospheric mass density ρ 0 at altitude z 0 , divided by the mass density ρ at altitude z, is taken instead of the dimensionless relative altitude Equation (34a), as the dependent variable
η ρ 0 ρ ( z ) = exp z z 0 = e ζ ;
this implies:
ζ = ( log η ) = η η ,
ζ = η η η η 2 ,
so that the exponential nonlinearity in the equation of motion (35) is replaced by algebraic nonlinearities
( 1 τ ) η η η η 2 + a η 2 b η 3 = f η 2 ,
viz. quadratic for the drag term on the right-hand side of Equation (41) and cubic for all other terms on the left-hand side of Equation (41), viz. inertia, gravity, and thrust. This will imply that drag effects will first appear at higher order than inertia, gravity and thrust.
The solution of the nonlinear differential equation (41) will be sought as a power series
η ( τ ) = n = 0 A n τ n + υ ,
with exponent υ and coefficients A n to be determined. Recalling Equations (39) and (34b), the trajectory is specified
exp z z 0 = n = 0 A n c ( t t 0 ) m 0 υ + n
by altitude (represented by the ratio of initial mass density to varying atmospheric mass density) as a function of time (represented by burned propellant mass relative to initial mass). The condition Equation (33b) of initial altitude z = z 0 at time t = t 0 , implies Equations (44a,b):
A 0 = 1 ,
υ = 0 ,
η ( τ ) = 1 + n = 1 A n τ n ,
specifying in Equation (44c) the exponent and first coefficient of Equation (42).
Using again Equations (34b) and (39), the trajectory Equation (44c) specifies:
exp z z 0 = 1 + n = 1 A n c ( t t 0 ) m 0 n ,
and by differentiation with regard to time, the velocity z ˙ in:
z ˙ exp z z 0 = c m 0 n = 1 n A n c ( t t 0 ) m 0 n 1 .
The initial condition Equation (33b) for the velocity z ˙ = v 0 at time t 0 specifies the coefficient A 1 :
v 0 = A 1 c m 0 = A 1 v * ,
v * c m 0 ,
as the ratio of initial velocity v 0 to the reference velocity v * , defined by Equation (47b). In order to interpret the latter note that the reference velocity v * corresponds Equation (48a) to a distance equal to the atmospheric scale height covered in a reference time t * :
v * = t * ,
t * m 0 c > t 1 t 0 ,
where the reference time is the initial rocket mass divided by the propellant burn rate; thus, the reference time would be the burn time if the rocket mass was entirely propellant, and must exceed the real burn time.

3.3. Leading Coefficients of Power Series Expansion

Using Equation (47a), the power series Equation (44c) for the trajectory becomes
η ( τ ) = 1 + v 0 v * τ + n = 2 A 0 τ n ,
and is considered next up to order N + 2 ; since the equation of motion (41) involves at most second-order derivatives, it becomes a polynomial of degree N in τ which is accurate only to order τ N and whose coefficients involve A 2 , A 3 , , A N + 2 ; equating to zero the coefficients of the N + 2 successive powers of τ , viz. τ 0 = 1 , τ , , τ N + 2 yields a set of N + 1 nonlinear algebraic equations, which can be solved for the N + 1 unknowns A 2 , , A N + 2 since A 0 , A 1 are already known from Equations (44a) and (47a). This process is illustrated by substituting the series solution for the trajectory Equation (49) for N = 2 , i.e., up to N + 2 = 4 the fourth-order in τ ,
η ( τ ) = 1 + v 0 τ v * + A 2 τ 2 + A 3 τ 3 + A 4 τ 4 + O ( τ 5 ) ,
into the Equation (41) up to order 2,
f v 0 v * + 2 A 2 τ + 3 A 3 τ 2 2 = b 1 + v 0 τ v * + A 2 τ 2 3 + ( 1 τ ) 1 + v 0 τ v * + A 2 τ 2 1 + v 0 τ v * + A 2 τ 2 2 A 2 + 6 A 3 τ + 12 A 4 τ 2 v 0 v * + 2 A 2 τ + 3 A 3 τ 2 2 + a 1 + v 0 τ v * + A 2 τ 2 2 + O ( τ 3 ) ,
which is valid only to O ( τ 3 ) because some cubic and higher-order terms were omitted.
Equating to zero the term of Equation (51) independent of τ specifies
2 A 2 = b a + ( 1 f ) ( v 0 / v * ) 2 ,
which is the coefficient of τ 2 in Equation (50) of the trajectory. Equating to zero the coefficient of τ in the equation of motion (51) leads to:
6 A 3 = ( 3 b 2 a ) v 0 v * + 2 A 2 ( 1 2 f ) v 0 v * + 1 v 0 v * 2 A 2 + a v 0 v * 2 ,
which on substitution of Equation (52) specifies the coefficient A 3 of τ 3 in the Equation (50) of the trajectory, viz.,
6 A 3 = b + ( b a ) ( 3 2 f ) v 0 v * f v 0 v * 2 + ( 1 2 f + 2 f 2 ) v 0 v * 3 .
Equating to zero the coefficient of τ 2 in the equation of motion (51) leads to
12 A 4 = 2 f 3 A 3 v 0 v * + 2 A 2 2 + 3 b A 2 + v 0 v * 2 3 a A 2 1 v 0 v * v 0 v * + 3 A 2 v 0 v * 2 + 6 A 3 1 v 0 v * v 0 v * 3 ,
which on substitution of A 2 , A 3 from Equations (52, 53a) specifies the coefficient A 4 of τ 4 in the equation of the trajectory (50), as a polynomial of degree four in v 0 / v * , with coefficients involving a , b , f . In a similar way, the vanishing of the coefficient of τ N in the equation of motion (41), specifies the coefficient A N + 2 of τ N + 2 in the equation of the trajectory (49); as can be seen from Equations (52, 53a, 54), the coefficient A N + 2 is a polynomial of v 0 / v * of degree N + 2 , whose coefficients involve the dimensionless gravity Equation (36a) = Equations (38c,d), thrust Equation (36b), and drag Equation (36c) = Equations (38a,b) parameters.

4. Alternative Methods for Short Times or Long Times

The dimensionless time variable is taken to be, instead of mass fraction of burned propellant (Section 3), which increases with time, the residual mass fraction (Section 4), which reduces with time (Section 4.1), and is minimum at burn-out, and thus is most convenient to calculate burn-out altitude and velocity (Section 5), once the series solution (Section 4.2) has been determined (Section 4.3).

4.1. Residual Mass Fractions as Time Variable

The dimensionless time variable is taken to be, instead of mass fraction of burned propellant Equation (34b), the residual mass fraction
1 μ m ( t ) / m 0 = 1 c ( t t 0 ) / m 0 m 1 / m 0 ,
which varies from unity μ = 1 at initial time t = t 0 to the ratio of empty mass at burn-out to total initial mass, when it is usually much smaller than unity. Using the residual mass fraction Equation (55) as an independent variable, and the same dependent variable Equation (34a), the equation of motion (32) becomes
μ ζ + a = b f ζ 2 e ζ ,
instead of Equation (35), where prime now denotes derivative with regard to μ , viz. ζ d ζ / d μ , and the dimensionless gravity Equation (36a) ≡ Equations (38c,d), thrust Equation (36b) and drag Equation (36c) ≡ Equations (38a,b) coefficients are the same.
Since the independent variable, viz. the residual mass fraction decreases with time, the dependent variable should decrease with altitude, e.g., the inverse of Equation (39) is used, viz. the ratio of atmospheric mass density ρ at altitude z to initial mass density ρ 0 at altitude z 0 is used:
ξ 1 / η = ρ ( z ) / ρ 0 = exp z z 0 = e ξ ,
as dimensionless altitude variable. Bearing in mind that:
ζ = ( log ξ ) = ξ ξ ,
ζ = ξ ξ + ξ ξ 2 ,
the exponential nonlinearity in the equation of motion (56) is replaced by algebraic nonlinearities
f ξ 2 ξ = b ξ 2 + μ ( a ξ 2 + ξ 2 ξ ξ ) ,
which are cubic for the drag term on the l.h.s. and quadratic for all the other terms on the right-hand side, viz. the thrust, gravity and inertia terms; this is the reverse of the situation before with Equation (41).
The solution of the equation of motion (59), i.e., the trajectory, is sought as the ratio of atmospheric mass densities Equation (52) expressed by a power series of residual mass fraction (55), viz.,
ξ = n = 0 B n μ n + ϑ
with exponent ϑ and coefficients B n to determined. The altitude is given explicitly Equation (57) as a function of time Equation (55) by
exp z z 0 = n = 0 B n 1 c ( t t 0 ) m 0 n + ϑ ,
and the velocity is obtained by differentiation with regard to time
z ˙ exp z z 0 = c m 0 n = 0 ( n + ϑ ) B n 1 c ( t t 0 ) m 0 n + ϑ 1 .
Substituting the initial conditions Equations (33a,b) respectively in Equations (61) and (62) leads to:
1 = n = 0 B n ,
v 0 = v * n = 0 ( n + ϑ ) B n ,
where the reference velocity Equation (47b) ≡ Equations (48a,b) appears again in Equation (63b); thus, the initial conditions Equations (63a,b) in method II involve all coefficients instead of just two in method I in Equations (44a) and (47a). In addition, they do not determine the exponent ϑ in Equation (60) instead of υ in Equation (44b).

4.2. Atmospheric Mass Density as a Function of Residual Mass

If the trajectory Equation (60) up to order N + 2 in μ is substituted in the equation of motion (59), since the latter involves at most second-order derivatives in ζ , it is accurate as a polynomial of degree N in μ . Equating to zero the coefficients of μ 0 1 , μ , , μ N leads to a set of nonlinear algebraic equations for the coefficients B 0 , B 1 , , B N ; the initial conditions Equations (63a,b) then specify B N + 1 and B N + 2 , completing the determination of the coefficients of the series up to order N + 2 . The exponent ϑ is determined at the outset by equating the lowest powers in the equation of motion. The process is illustrated next for N = 2 , thus specifying for method II the exponent and coefficients of the series up to order N + 2 = 4 , viz. B 0 , B 1 , B 2 , B 3 , B 4 ; this corresponds to the determination in method I of A 0 , A 1 , A 2 , A 3 , A 4 in Equations (44,47a,52,53b,54), but leads to simpler algebra in the present case, as shown next.
In order to determine the exponent ϑ , only the leading power in the trajectory Equation (60), viz. Equation (64a):
ξ = B 0 μ ϑ 1 + O ( μ ) ,
f B 0 ϑ 2 μ 3 ϑ 2 = b μ 2 ϑ + a μ 2 ϑ + 1 + ϑ 2 μ 2 ϑ 1 ϑ ( ϑ 1 ) μ 2 ϑ 1 ,
is substituted in the equation of motion (59), leading to Equation (64b) where only the lowest powers are valid, viz:
B 0 f ϑ μ 3 ϑ 2 = μ 2 ϑ 1 .
It follows from Equation (65) that the exponent must satisfy 3 ϑ 2 = 2 ϑ 1 , leading to Equation (66a), and hence the coefficient Equation (66b):
ϑ = 1 ,
B 0 = 1 / f .
Using the coefficients C n defined by Equation (66c) instead of the B n ,
C n B n f ,
leads to
f ξ = μ + n = 1 C n μ n + 1 ,
as Equation (60) of the trajectory.
In the first method I of trajectory calculation (Section 3), the power series solution Equation (49), was calculated explicitly Equation (50) up to the fifth term inclusive, and the same will be done with the second method II of trajectory calculation (Section 4), for the power series solution Equation (67), viz.:
f ξ ( μ ) = μ + C 1 μ 2 + C 2 μ 3 + O ( μ 4 ) ,
which is taken only to O ( μ 3 ) because, once C 1 , C 2 are known, then the initial conditions Equations (63a,b) with Equation (66a,b,c):
f = 1 + C 1 + C 2 + C 3 + C 4 ,
f v 0 v * = 1 + 2 C 1 + 3 C 2 + 4 C 3 + 5 C 4 ,
specify C 3 , C 4 . The equation of trajectory (68) is substituted to the same order O ( μ 3 ) in the equation of motion (59), viz.
μ ( 1 + C 1 μ + C 2 μ 2 ) ( 1 + 2 C 1 μ + 3 C 2 μ 2 ) 2 = b μ 2 ( 1 + C 1 μ ) 2 + μ a μ 2 + ( 1 + 2 C 1 μ + 3 C 2 μ 2 ) 2 μ ( 1 + C 1 μ ) ( 2 C 1 + 6 C 2 μ ) ,
where Equation (70) is valid only to order μ 3 , since some terms of order μ 4 were already omitted.

4.3. Determination of Exponent and Coefficients of the Series Solution

Note that, in the second method II, the exponent and first coefficient Equations (66a,b) are determined at the leading order Equation (65); in addition, the initial conditions Equations (63a,b) supply two extra coefficients ( B N 1 , B N ) or ( C N 1 , C N ) by Equation (66c), so that the calculation of the series solution Equation (60) by method II requires only an expansion to order μ N 2 in the equation of motion (59) for an accuracy of order μ N . In the first method I, the calculation of the power series Equation (49) to order μ N requires calculating ( N + 1 ) coefficients in the equation of motion (41), which must be taken to order μ N ; this is two orders higher than μ N 2 in method II, so that the latter involves much less tedious algebra for the same accuracy. This can be confirmed for N = 4 , comparing Equation (51), which determines the coefficients ( A 2 , A 3 , A 4 ) for method I, with Equation (70), which determines the coefficients ( C 1 , C 2 ) for method II, with ( C 3 , C 4 ) specified by Equations (69a,b). The resulting equations for the determination of ( A 2 , A 3 , A 4 ) in method I are more complex Equations (52, 53b, 54) than those for method II, as will be shown next, confirming that the latter is algebraically simpler to implement.
The coefficient of the lowest power μ in Equation (70) leads to an identity 1 = 1 because it was already used Equation (65) to determine the exponent Equation (66a) and first coefficient Equation (66b) of the series Equation (68). Equating to zero the remaining coefficients of μ 2 and μ 3 in Equation (70), leads respectively to:
3 C 1 = b ,
7 C 2 = a 2 b C 1 6 C 1 2 = a ,
where Equation (71a) determines explicitly the coefficient C 1 , which may be substituted in Equation (71b) to specify explicitly C 2 . Using Equations (71a,b) in the initial conditions Equations (69a,b) leads to the system of equations:
C 3 + C 4 = 1 f a 7 + b 3 ,
4 C 3 + 5 C 4 = 1 f v 0 v * 3 a 7 + 2 b 3 ,
which specify the remaining two coefficients:
C 3 = 4 f 5 v 0 v * 2 a 7 + b ,
C 4 = 3 + f 4 v 0 v * + a 7 2 b 3 .
The coefficients Equations (71a,b; 73a,b) determine the series solution for the trajectory Equation (67) with five terms up to order μ 4 .
Only the first three terms Equation (68) will be written explicitly because, by Equations (69a,b), they already involve:
f ξ ( μ ) = μ b μ 2 + ( a 4 b 2 ) μ 3 ,
all three dimensionless parameters of the trajectory, viz. the gravity Equation (36a), thrust Equation (36b) and drag Equation (36c) factors. Using also Equations (55) and (57), the first three terms of the equation of the trajectory are written explicitly,
C D 2 exp z z 0 = m 0 ρ 0 S 1 c ( t t 0 ) m 0 + m 0 2 T ρ 0 S c 2 2 1 c ( t t 0 ) m 0 2 m 0 3 g ρ 0 S c 2 2 4 m 0 3 T 2 ρ 0 S c 4 3 1 c ( t t 0 ) m 0 3 + ,
which can be simplified in terms of the variable mass Equation (2a) and reference mass Equation (38a), viz.,
C D 2 exp z z 0 = m ( t ) m * + T c 2 [ m ( t ) ] 2 m * g c 2 4 T 2 c 4 2 [ m ( t ) ] 3 m * + ,
in order for Equations (75, 76) to meet the initial conditions Equations (33a,b), the next two terms C 3 μ 3 + C 4 μ 4 involving Equations (73a,b) would have to be included. The result Equation (76) will be analyzed in more detail subsequently (Section 6.2); here, it is noted that the left-hand side of Equation (76) must be positive, and the right-hand side shows that the only positive term, which contributes to the climb, is thrust, whereas gravity enters negatively, opposing the climb.
Methods I and II of calculation of rocket trajectories share the following assumptions: (i) vertical climb; (ii) constant propellant flow rate; (iii–iv) constant thrust aligned with the rocket axis; (v) zero angle of attack, hence no lift; (vi–vii) drag a quadratic function of the velocity in the absence of wind; (viii) mass density decaying exponentially with altitude for an isothermal atmospheric layer; (ix) uniform gravity over an altitude range small compared to radius of the Earth; (x–xi) negligible Coriolis and centrifugal forces due to the rotation of the Earth; and (xii) single stage rocket. Some of these restrictions can be easily removed, e.g., for a multi-stage rocket using the final condition of one stage as the starting condition for the next stage. The two methods differ in that (i) method I uses the dimensionless propellant mass consumption divided by the initial mass as variable, which is smaller for a short time after lift-off so that the series solution is more accurate with less terms; (ii) method II uses the dimensionless residual mass divided by the initial mass as variable, which is smaller closer to burn-out, so that again the series solution is more accurate with less terms. The short time after lift-off of method I and the short time before burn-out method II can be matched in the intermediate section of the trajectory. An alternative to method I for a short time after lift-off is method III using a power series of time. A complement to method II after burn-out is method IV for the unpowered ballistic trajectory. In the case of method IV, all terms of the series solution are obtained explicitly, whereas, for methods I, II, and III, the computation is illustrated for the first four to five terms.

5. Three Distinct Methods of Trajectory Calculation

The two preceding methods of trajectory calculation are compared for accuracy (Section 5.1), showing that method II is more suited for a long time (Section 4), e.g., burn-out, and method I for a short time (Section 3), e.g., lift-off, although both cover the full time range. Since method I has heavier algebra, a simpler third method III is developed (Section 5.2) also suited for a short time; a particular case with exact analytical solution as a power series with terms of all orders calculated explicitly, is the case of ballistic flight without thrust, after burn-out (Section 5.3).

5.1. Comparison of Accuracy for Short Times and Long Times

The trajectory Equation (45) calculated by (Section 3) altitude as a function of variable mass Equation (2a), viz.
exp z z 0 = 1 + n = 1 A n [ 1 m ( t ) / m 0 ] n ,
is a power series of burned propellant mass ratio, implying that it converges more quickly in a short time, when the mass is still relatively close to lift-off mass Equation (78a):
1 m ( t ) / m 0 ϵ ,
E N n = N | A n | ϵ n ,
so that the error after N terms has an upper bound Equation (78b). If a sufficiently large number of terms Equation (79a) of the converging series Equation (77) is calculated, it can be expected that the coefficients decay at a rate Equation (79b):
n N ,
| A n + 1 / A n | r < 1 ,
| A n + N | < r n | A n | ,
implying Equation (79c) so that an upper bound for the error in Equation (78b) is
E n | A N | ϵ N n = 0 r n ϵ n = | A N | ϵ N 1 r ϵ .
Thus, the accuracy is high for short time ϵ 1 , and increases rapidly as ϵ N when adding more terms to the series.
The worst case for method I of calculation of the trajectory Equation (77) is the longest time, viz. burn-out, which occurs at altitude z 1 specified by
exp z 1 z 0 = 1 + n = 1 A n ( 1 m 1 / m 0 ) n ,
where, if the residual mass fraction is small m 1 m 0 then ( 1 m 1 / m 0 ) n 1 , corresponding to ϵ = 1 in Equation (80), so that convergence depends on the coefficients Equation (79b) alone. This was the reason for the consideration of method II which (Section 4) specifies the trajectory by Equation (61) where the variable mass Equation (2a) can also be substituted,
C D 2 m 0 m * exp z z 0 = 1 + n = 0 C n m ( t ) m 0 n ,
where Equations (66a,b,c) and (38b) were used. It follows that Equation (82) converges best for long times, when the variable mass is closest to residual mass, and is a smaller fraction of initial mass.
Thus, Equation (82) is most suited to calculate burn-out conditions,
C D 2 m 0 m * exp z 1 z 0 = n = 0 C n m 1 m 0 n ,
with an upper bound for the error after N terms,
E ¯ N n = N | C n | m 1 m 0 n = m 1 m 0 N n = 0 | C N + n | m 1 m 0 n .
Assuming that, for sufficiently high-order Equation (85a), the coefficients have a decay rate Equation (85b):
n N ,
| C n + 1 / C n | r ¯ < 1 ,
| C N + n | r ¯ n | C n | ,
then Equation (85c) implies that an upper bound for the error Equation (84) is
E ¯ N m 1 m 0 N | C N | 1 r ¯ ( m 1 / m 0 ) ;
for small residual mass fraction m 1 m 0 , the term in the denominator can be omitted, and accuracy increases rapidly with order N. Replacing in Equation (86) m 1 by m ( t ) , it is clear that the accuracy is worst for short time m m 0 when Equation (86) reduces to | C N | / ( 1 r ¯ ) , showing that convergence depends on the coefficients alone. It has been shown that methods I (Section 3) and II (Section 4) of trajectory calculation are complementary, in that one works best when the other is worst.

5.2. Direct Calculation of the Taylor Series Expansion

A significant difference between the two preceding methods is that: (a) method II which applies best to long times has moderate algebra (Section 4); (b) method I has heavier algebra (Section 3), but is preferable from the point of view of accuracy for short-time. This suggests the introduction as an alternative to method I, of a method III, which converges rapidly for short times, but has lighter algebra. Note that since methods I and III converge best for short times, say t 0 < t < t 2 t 1 , the parameters of the rocket and the atmosphere can be easily updated at each successive time interval, e.g., t 2 < t < 2 t 2 or in general N t 2 < t < ( N + 1 ) t 2 < t 1 ; method II converges best close to burn-out t 0 < t < 2 t 2 but assumes that the same rocket and atmospheric parameters apply to the whole flight. Since method III aims to converge most rapidly for short times, a Taylor series expansion of altitude as a function of time may be used:
z ( t ) = n = 0 D n ( t t 0 ) n ,
n ! D n = z ( n ) ( t 0 ) ,
where the initial conditions Equations (33a,b) specify the first two coefficients
z ( t ) = z 0 + v 0 ( t t 0 ) + n = 2 D n ( t t 0 ) n ,
and the remaining coefficients can be calculated directly from the equation of motion (32).
The equation of motion (32) specifies immediately the initial acceleration, using the reference mass Equation (38b),
z ¨ 0 z ¨ ( t 0 ) = g + T m 0 C D 2 m * m 0 v 0 2 ,
which determines the coefficient of the third term of the trajectory Equation (88). Differentiating the equation of motion (32) once with regard to time,
[ m 0 c ( t t 0 ) ] z c ( z ¨ + g ) = C D S ρ 0 2 z ˙ 3 / 2 z ˙ z ¨ exp z z 0 ,
and choosing the initial time
m 0 z 0 = c ( z ¨ 0 + g ) + C D 2 m * ( v 0 / ) ( v 0 2 / 2 z ¨ 0 )
specifies the fourth coefficient
z 0 = c T m 0 2 + C D m * m 0 v 0 g T m 0 v 0 c 2 m 0 + v 0 2 1 + C D 2 m * m 0 ,
where Equation (89) was used. All the remaining coefficients Equation (87b) of the trajectory Equation (87a) can be calculated the same way.
In the case of lift-off with zero initial velocity at zero altitude, the coefficients Equations (89) and (90c) simplify respectively to:
z 0 = 0 = v 0 : z ¨ 0 = g + T m 0 ,
z 0 c T m 0 2 ;
substituting in the trajectory Equations (87a,b) up to fourth-order:
z ( t ) = z 0 + v 0 ( t t 0 ) + ( z ¨ 0 / 2 ) ( t t 0 ) 2 + ( z 0 / 6 ) ( t t 0 ) 3 ,
and putting initial time, altitude and velocity all equal to zero leads to
t 0 = 0 : z ( t ) = t 2 2 g + T m 0 1 + c t 3 m 0 ,
showing that the initial acceleration T / m 0 g is corrected to next order t 3 by the propellant flow rate, and drag would appear only at higher order t 4 . The three preceding methods I, II and III of trajectory calculation using series expansions can be replaced by an analytic closed form solution in the case of free flight with zero thrust after burn-out; the latter solution is known in the case of constant atmospheric density [3], and is extended next to an isothermal atmosphere, in a form which is suited as a sequel to method II (or I or III) after burn-out.

5.3. Peak Altitude in Post Burnout Flight

The assumption of zero initial altitude and velocity used in Equation (93) applies only to the first stage of a multi-stage rocket; it could also be used to simplify the preceding results for methods I to III. The latter were presented for any initial velocity and altitude in order to apply to all stages of multistage rocket. Thus, any of the methods I to III can be applied to calculate burnout altitude z 1 and velocity v 1 z ˙ 1 , although: (i) method II is the most accurate for constant rocket and atmospheric parameters; (ii) methods I and III may be preferable otherwise. After burnout, there is no thrust or propellant consumption, and thus the rocket mass is constant and equal to the residual mass m 1 , simplifying the equation of motion (32) to
m 1 ( z ¨ + g ) = 1 2 C D S ρ 1 z ˙ 2 exp z z 1 ,
where:
m 1 = m 0 c ( t 1 t 0 ) ,
ρ 1 ρ 0 exp z 1 z 0 ,
the atmospheric mass density at burnout is (95b). Using the identity
z ¨ d z ˙ d t = d z ˙ d z d z d t = z ˙ d z ˙ d z ,
the equation of motion (94) can be written as
d d z z ˙ 2 + 2 g = C D ρ 1 S m 1 z ˙ 2 exp z z 1 ,
in terms of z ˙ 2 , which is twice the kinetic energy per unit mass.
Using as dependent variable Equation (98a), a dimensionless form of the kinetic energy per unit mass:
Ψ z ˙ 2 2 g ,
φ exp z z 1 = ρ ( z ) ρ 1 ,
and using as dimensionless altitude variable Equation (57) the ratio of mass densities Equation (98b) at burnout altitude z 1 , and at altitude z, the equation of motion (97) becomes
C D ρ 1 S m 1 φ Ψ + 1 = d Ψ d z = d Ψ d φ d φ d z = φ d Ψ d φ ,
which involves Equation (100a):
d Ψ d φ 2 h Ψ = 1 φ ,
h = C D ρ 1 S 2 m 1 ,
a single dimensionless parameter Equation (100b), viz. the drag parameter Equation (36c), calculated for the residual rocket mass Equation (95a) and atmospheric mass density Equation (95b) at burn-out altitude. Note that the dimensionless dependent variable Ψ = H / in Equation (98a), is the ratio by the scale height of H = Ψ = z ˙ 2 / ( 2 g ) , which is the height reached by a body launched vertically at initial velocity z ˙ against a uniform gravity field g. The solution of the equation of motion Equation (100a) is sought in the form of Equation (101a):
Ψ ( φ ) = e 2 h φ Φ ( φ ) ,
d Φ d φ = φ 1 e 2 h φ ,
where substitution of Equation (101a) in Equation (100a) shows that the function Φ satisfies Equation (101b). Expanding the exponential in Equation (101b) in power series:
d Φ d φ = 1 φ + n = 1 ( 2 h ) n n ! φ n 1 ,
and integrating:
Φ ( φ ) = Φ 1 + log φ + n = 1 ( 2 h ) n n ! n φ n ,
specifies the function Φ ( φ ) to within an arbitrary constant Φ 1 .
Noting that Equation (98b) and Equation (100b) imply
exp ( 2 h φ ) = exp C D S m 1 ρ 1 exp z z 1 = exp C D S m 1 ρ ( z ) ,
and substituting together with Equations (102b, 98b) in Equations (101a, 98a) specifies the ascent velocity z ˙ at arbitrary altitude z,
z ˙ 2 2 g = exp C D S m 1 ρ ( z ) Φ 1 z z 1 + n = 1 C D S m 1 ρ ( z ) n 1 n ! n ,
where the constant Φ 1 can be determined from the initial ascent velocity z ˙ 1 at burn-out altitude z 1 , viz.
Φ 1 = z ˙ 1 2 2 g exp C D S ρ 1 m 1 n = 1 C D S ρ 1 m 1 n 1 n ! n .
Substituting Equation (105) in Equation (104) specifies the ascent velocity as
z ˙ 2 exp C D S m 1 ρ ( z ) z ˙ 1 2 exp C D S m 1 ρ 1 = 2 g ( z z 1 ) + 2 g n = 1 C D S m 1 n 1 n ! n [ ρ ( z ) ] n ρ 1 n ,
which becomes an identity 0 = 0 for z = z 1 , and shows that the ascent velocity z ˙ decreases with z, because the right-hand side of Equation (106) is negative. The ascent velocity vanishes z ˙ 2 = 0 at the peak altitude z 2 which is a root of
z ˙ 2 = 0 : z 2 z 1 = z ˙ 1 2 2 g exp C D S ρ 1 m 1 + n = 1 C D S ρ 1 m 1 n 1 n ! n exp z 2 z 1 1 n .
Note that the first term on the right-hand side of Equation (107) is positive, and also all the terms of the series are positive because: (i) for even order n = 2 , 4 , 5 , both powers are positive; (ii) for n odd n = 1 , 3 , 5 , the expressions in square and curved brackets are negative, so the product is still positive. The series in Equation (107) decays rapidly, and a few terms specify a transcendental equation whose root z 2 is the peak altitude.

6. Computation of a Powered Followed by a Ballistic Ascent

A numerical example of the preceding methods of trajectory calculation is given (Section 6.1) by considering: (i) a powered first stage, for which any of the methods I to III can be used to specify the burn-out altitude z 1 and burn-out velocity z ˙ 1 in (Section 6.2); (ii) these are the initial conditions for a ballistic ascent (Section 6.3) for which method IV specifies the peak altitude z 2 .

6.1. Input Data for Trajectory Calculation

The equations describing the powered trajectory of the rocket are Equation (45) for the method I, Equation (61) for the method II and Equation (88) for the method III. All three equations have a similar form: (i) on the right-hand side a function of time involving the rocket parameters in the coefficients; (ii) on the left-hand side a function of altitude specified by the atmospheric density profile. Thus the methods I, II and III apply to any atmospheric density profile. For example, in the method I using Equations (34b), (39) and (44b) in Equation (42) the trajectory is given by
ρ 0 ρ ( z ) = 1 + i = 1 A i c ( t t 0 ) m 0 i ,
for any atmospheric density profile. In the particular case of an isothermal atmosphere Equation (39) substituted in Equation (108) leads to Equation (45). For a polytropic atmosphere [13] with a constant temperature gradient Equation (109a) where θ is the temperature, θ 0 the temperature at level z 0 = 0 and 0 the lengthscale
d θ d z = θ 0 1 :
ρ ( z ) ρ 0 = 1 z 1 1 / ( n 1 ) ,
the profile of mass density is Equation (109b), where n is the polytropic index. Thus the trajectory of the same rocket in a polytropic atmosphere is given substituting Equation (109b) in Equation (108) leading to Equation (110),
1 z 0 1 / ( n 1 ) = 1 + i = 1 A i c ( t t 0 ) m 0 i ,
instead of Equation (45) for an isothermal atmosphere.
The lowest layer of the atmosphere of the Earth is polytropic with the temperature reducing linearly from T 0 = 288   K at zero altitude through the troposphere to T 2 = 216.5   K at the tropopause z 2 = 11   k m . The temperature then remains constant in the stratosphere that is isothermal up to the altitude z 3 = 25   k m . The preceding data corresponds to the International Standard Atmosphere (ISA) and can be shifted in temperature, e.g., ISA + 10   ° C . Other atmospheric models may be used extending to higher altitudes. For a detailed trajectory calculation the trajectory equation should be applied for each layer, e.g., in model I starting with Equation (110) for the troposphere 0 z 11   k m , followed by Equation (45) for the stratosphere 11   k m z 25   k m , and so on. Since the purpose of the present Section is to illustrate the effect of rocket parameters on the trajectory a simplified “average” atmospheric model is used. The scale height for an isothermal atmosphere is given by = c 0 2 / ( γ g ) where the sound speed is c 0 = 340   m / second at sea level and the adiabatic exponent is γ = 1.4 for a diatomic perfect gas like air that is made mostly of nitrogen N 2 and oxygen O 2 . Using g = 9.81   m / s 2 leads to a scale height 0 = 8.4   k m . The constant temperature gradient between sea level θ 0 = 288   K and θ = 231.5   K at the tropopause z 2 = 11   k m leads by Equation (109a) to a lengthscale
1 = θ 0 d θ / d z = θ 0 ( θ 2 θ 0 ) / z 2 = z 2 1 θ 2 / θ 0 ,
equal to 1 = 44.3   k m . The arithmetic average of the scale height and lengthscale = ( 0 + 1 ) / 2 = 26   k m is taken as the scale height for a simplified isothermal atmosphere over the altitude range of rocket flight.
The preceding methods I to IV of trajectory calculation can be applied numerically, e.g., in small steps in a computer code, for maximum accuracy. For the purpose of illustrating the main features of each method, it is sufficient to use only a few discretization points, allowing a calculation by hand, which makes the main points more obvious. The trajectory could be illustrated alternatively and equivalently by plots or tables. Since tables are used to indicate the rocket data, they are also adopted to indicate the trajectory. The rocket data concerns the Ariane 5 that is currently the most reliable and widely used heavy launcher. The data in the Table 2 concerns the last model of the family, namely the Ariane 5 ECA [14] for which a typical booster release altitude can be found in [15]. The rocket consists of two stages and two boosters with a total lift-off mass m 0 = 777   t , where 1 t = 1000   k g ; the powered flight is considered from lift-off t 0 = 0   s to the time t 1 = 140   s of booster separation, and it is necessary to calculate the total mass m 1 at this time. Each booster has a total mass m 2 = 273   t and empty mass m 3 = 33   t , corresponding to a propellant load m 4 = m 2 m 3 = 240   t . This propellant is spent by the time t 1 = 140   s for both boosters reducing the total rocket mass by m 5 = 2 m 4 = 480   t . The core first stage has a gross mass m 6 = 184.7   t and empty mass m 7 = 14.7   t , and thus carries a mass m 8 = m 6 m 7 = 170   t of propellants and burns for t 2 = 540   s . Assuming a constant flow rate, the mass of propellants consumed up to the booster separation time t 1 = 140   s is m 9 = m 8 t 1 / t 2 = 0.259 m 8 = 44.07   t . The total propellant consumption of the two boosters m 4 = 240   t and of the core first stage m 9 = 44.07   t adds to m 10 = m 4 + m 9 = 284.07   t over a time t 1 = 140   s leading to an average fuel consumption c = m 10 / t 1 = 2.029   t / s . The lift-off of mass m 0 = 777   t is reduced by the propellant consumption m 10 = 284.07   t to m 1 = m 0 m 10 = 492.9   t at booster separation time t 1 = 140   s . The thrust of the core stage T 1 = 1390   k N and of the two boosters T 2 = 7080   k N add to a total vacuum thrust of T = T 1 + 2 T 2 = 15550   k N that is adopted as a constant value during the first phase of boosted flight. The core stage has a radius R 1 = 2.7   m and cross-sectional area S 1 = π R 1 2 = 22.9   m 2 and each booster has a radius R 2 = 1.53   m and cross-sectional area S 2 = π R 2 2 = 7.354   m 2 for a total cross-sectional area S = S 1 + 2 S 2 = 37.61   m 2 . The input data (Table 2) consists of three sets: (i) ‘environmental’ data, viz. the acceleration of gravity, the sea level mass density and atmospheric scale height; (ii) ‘propulsion’ data, viz. thrust, propellant flow rate and mass, and burn-time; (iii) ‘aerodynamic’ data, viz. lift-off mass, cross-sectional area and drag coefficient. Some of the data, like the aerodynamic coefficients, should be updated along the flight path for acceptable accuracy; since the aim here is illustration only, a constant data is used to calculate: (iv) the reference mass Equation (38a); (v) the weight Equation (36a), thrust Equation (36b), and drag Equation (38b) parameters; (vi–vii) the reference time Equation (48b) and velocity Equation (47b). The six reference parameters at the bottom of the Table 2 specify the performance of a rocket; two different rockets with the same six parameters would fly the same trajectory.

6.2. Altitude and Ascent Velocity in Powered Flight

The data in the Table 2 can be applied to any of the three methods I, II and III of calculation of a powered trajectory. In all methods for the purpose of illustration, it is sufficient to use a four-point discretization (Table 3), i.e., calculate the altitude z and vertical velocity z ˙ at four equally spaced fractions of the booster burn-time t = t 1 / 4 , t 1 / 2 , 3 t 1 / 4 , t 1 viz. t = 35   s , 70   s , 105   s , 140   s . The same data is used all points, viz. the coefficients of the corresponding series expansions in powers of time. For method I, the altitude is specified as a function of time by Equation (45), viz.,
I : exp [ z ( t ) / ] = 1 + n = 1 A n ( t / t * ) n ,
where t 0 = 0 and z 0 = 0 at lift-off, and the reference time was used Equation (48b). The first five coefficients of the series expansion are given by Equation (44a) for A 0 , Equation (47a) for A 1 , Equation (52) for A 2 , Equation (53b) for A 3 and Equation (54) for A 4 ; these involve the mass Equation (36a), thrust Equation (36b) and drag Equation (38b) parameters and reference velocity Equation (47b) and time Equation (48b) that are already available (Table 2) and are used in the trajectory calculation in Table 3. The expressions are simplified by the zero initial velocity v 0 = 0 of the first stage:
{ A 0 , A 1 , A 2 , A 3 , A 4 } = 1 , 0 , b a 2 , b 6 , b 12 + ( b a ) 2 24 ( 3 2 f ) ,
where the last coefficient was calculated from Equation (54):
12 A 4 = 4 f A 2 2 + 3 ( b a ) A 2 + 6 A 3 ,
using the third and fourth components of Equation (113). Substitution of Equation (113) in Equation (112) leads to
exp ( z / ) = 1 + b a 2 t t * 2 + b 6 t t * 3 + A 4 t t * 4 ,
or equivalently using the values in the Table 2,
exp ( z / ) = 1 + 28.77 t t * 2 + 18.81 t t * 3 + 391.4 t t * 4 .
For a time t 1 = 140   s or t 1 / t * = 0.3656 of booster separation the altitude exp ( z 1 / ) = 12.76 is z 1 = 2.546 = 66.2   k m in agreement with the published [15] flight data.
Table 3 also includes the vertical velocity, calculated in method I by Equation (46) with Equation (47b), viz.,
I : z ˙ ( t ) = v * e z / n = 1 n A n ( t / t * ) n 1 ,
under the same assumptions as Equation (112). Using the preceding coefficients in Equation (113) leads to the vertical velocity to the fourth order
v ( t ) = v * exp ( z / ) t t * b a + b 2 t t * + 4 A 4 t t * 2 ,
or equivalently
v ( t ) = v * exp ( z / ) t t * 57.54 + 56.43 t t * + 1565.5 t t * 2 ,
The exponential factor in Equation (119) has a strong effect in reducing the ascent velocity for altitudes in excess of the scale height z > , resulting in a rather low vertical burn-out velocity, because a vertical climb is inefficient in wasting energy against gravity. Following the vertical launch of a rocket, a gravity turn is used to tilt the velocity vector progressively away from the vertical, so that the vertical component of the velocity becomes a smaller fraction of the total velocity, and vanishes on insertion on a circular orbit. Thus, the vertical velocity component has a maximum between the zeros at launch and at orbital insertion. The Table 3 suggests that the peak of the vertical velocity component is about v 3 = 0.62   k m / s at altitude z 3 = 45.5   k m at a time t 3 = 105   s after launch. These calculations assume a vertical trajectory rather than a gravity turn; the present methods I, II and III can be extended from one to two dimensions to model a gravity turn. The comparison of method I with methods II, III and IV proceeds for a vertical trajectory.
The calculation of altitude and ascent velocity as a function of time is compared next from the preceding method I (in Section 3), to method II (in Section 4) and method III (in Section 5.2). In method II, the altitude is specified by Equation (67) with Equation (57), viz.
II : f exp z ( t ) / = ( 1 t / t * ) 1 + n = 1 C n ( 1 t / t * ) n ,
where were used Equations (55) and (48b):
μ = 1 c t m 0 = 1 t t * .
The ascent velocity is specified by differentiation of Equation (120),
II : f z ˙ ( t ) = v * e z / 1 + n = 1 ( n + 1 ) C n ( 1 t / t * ) n .
Note that the power series of time in Equation (120) must add to a value in modulus less than | f | because e z / < 1 in method II, in contrast to a value larger than unity e z / > 1 in Equation (112) in method I. Thus, method I can be used for large time steps as seen before (Section 6.1, from lift-off to burn-out 0 t t 1 ; method II can only be used for small time steps close to burn-out 1 ϵ 1 t / t 1 1 , and thus should not be applied to the coarse time grid 1 t / t 1 = 1 / 4 , 1 / 2 , 3 / 4 , 1 in Table 3.
The method III specifies the altitude by Equation (88) and ascent velocity respectively by:
III : z ( t ) = n = 0 D n t n ,
z ˙ ( t ) = n = 2 n D n t n 1 ,
which have been simplified using z 0 = 0 = v 0 . This cancels the first two coefficients Equation (87b), and the next two are specified respectively by Equations (89) and (90c), simplified to Equations (91a,b), viz.:
2 D 2 = g + T / m 0 ,
6 D 3 = c T / m 0 2 .
The first non-zero coefficient Equation (124a) is the excess of thrust to weight ratio at lift-off over the acceleration of gravity, and the second Equation (124b) involves thrust again, but not drag:
z ( t ) = t 2 ( D 2 + D 3 t ) ,
z ˙ ( t ) = t ( 2 D 2 + 3 D 3 t ) .
The results in Table 3 show a rapid rise in altitude at later stages of flight and an increasing ascent velocity leading to a high value at burn-out, in contrast with method I. In method III, an exagerated burn-out altitude and velocity results from using only two terms of the series Equations (123a,b) whose coefficients do not yet include drag. It is clear that at least one more term in the series would be needed to include drag for more accurate results.

6.3. Peak Altitude in Unpowered Ballistic Flight

The calculation of the trajectory post burn-out can be made using method IV for which a complete series expansion with all coefficients explicit, Equation (106), is available (Section 5.3). If the purpose is to calculate peak altitude Equation (107), then only a few parameters are needed. The atmospheric mass density Equation (95b) at burn-out Equation (126a):
ρ 1 = ρ 0 e z 1 / ,
m * * = ρ 1 S ,
appears in the reference mass Equation (126b), which is similar to Equation (38a) calculated at burn-out altitude instead of lift-off. It appears in the drag parameter Equation (100b) similar to Equation (38b), but using the initial mass of the upper stage in Equation (127a):
ϑ C D m * * / m 1 2 h ;
E ¯ 0 z ˙ 1 2 / ( 2 g ) ,
in addition to Equation (126a) the remaining dimensionless parameter is the kinetic factor Equation (127b), and both appear together with the preceding in Table 4.
The parameters Equations (127a,b) also appear in the combination Equations (128a,b):
E 0 = E ¯ 0 exp ( ϑ ) ,
X ( z 2 z 1 ) / ,
when calculating the relative burn-out altitude Equation (128b), which is given by the exact formula (107), viz.,
X = E 0 + n = 1 ( ϑ ) n n ! n e X 1 n ;
this is a transcendental equation involving the parameters E 0 , ϑ , whose positive real roots X > 0 specify the burn-out altitude. In the present case, the parameter ϑ 1 is rather small as seen in Table 4, where the burn-out data from methods I and III were used to calculate the initial conditions for ballistic flight. The approximation of small ϑ simplifies Equation (129) to
X E 0 = ϑ e X 1 + ϑ 2 4 e X 1 2 + O ( ϑ 3 ) ,
to lowest order
ϑ 3 , ϑ X 2 1 : X = E 0 1 ϑ .
As an example is simulated an Ariane 5 ECA abort at the time t 1 = 140   s of booster separation; it is assumed that the boosters remain attached so the mass is m 1 = 500.4   t . The booster propellant is exhausted and if the core engine is cut-off there is no thrust leading to a ballistic flight with constant mass. The results in the Table 4 show that since method I yields a much lower burn-out velocity than method II, the altitude gain is much smaller, and the peak altitude lower, confirming the obvious fact that ballistic flight is very sensitive to initial conditions.
The preceding account was not intended to be an accurate trajectory calculation, as shown by the very coarse time grid in one-fourth of the booster burn time. This was a deliberate choice to test and probe in a hand calculation the points which need attention in an accurate computation, viz.: (i) for method I in Section 3, short time steps with updated data, more so for the ascent velocity Equation (117) than the altitude Equation (112); (ii) method II in Section 4 holds only if the sum of the series Equation (120) is less than unity, for very short time steps taken backwards in time; (iii) method III in Section 5.2 needs the inclusion of drag terms starting with D 5 in Equations (123a,b) to be accurate at longer times; and (iv) using accurate burn-out data, method IV in Section 5.3 can lead to an precise estimation of peak altitude in ballistic flight as a root Equation (128b) of the transcendental equation Equation (129), whose solution is simple Equation (130b) for small ϑ and moderate X. In spite of its simplicity, method I was able to estimate accurately the booster separation altitude of the Ariane 5 ECA as observed in flight. The lessons learned on the preceding methods should remain valid in their application in a subsequent work to the far more relevant case of a two-dimensional gravity turn; in connection with the latter, the two-point boundary-value problem (TPBVP) of given a lift-off condition meeting a desired burn-out condition can also be addressed; the TPBVP is not always feasible (e.g., if the rocket has insufficient performance), and its discussion is deferred to its proper context of a gravity turn in future work.

7. Conclusions

The four methods I to IV of calculation of rocket trajectories have been applied to a single isothermal atmospheric layer, and use “analytical” solutions of the equations of motion; their accuracy would be improved by dividing the atmosphere into layers as is currently done in the “numerical” solution of the equations of motion. The main difference between the current “numerical” computation of rocket trajectories and the ensemble of the four “analytical” methods introduced in the present paper is that, for the same number N of atmospheric layers, the “analytical” approach should be more accurate than the “numerical” approach because: (i) physically it takes into account the variation of atmospheric density in the solution of the equations of motion, instead of taking the atmospheric density as a constant in each layer and discretizing the differential equation; and (ii) mathematically, for each atmospheric layer or step, the “analytical” solution can be more accurate than a numerical approximation to the solution of the differential equation.
Thus, in the comparison between the “analytical” and the “numerical” approaches: (i) the “analytical” approach is more accurate for the division of the atmosphere in the same number of steps; (ii) for a given accuracy, the “analytical” approach will require a smaller number of steps or atmospheric layers for the same altitude range; (iii) an intermediate choice between (i) and (ii) is to use the “analytical” approach with a smaller number of steps than the “numerical” approach and still obtain higher accuracy. The current approach to the calculation of rocket trajectories can be called “numerical” because it is based on the discretization of the equations of motion and division of the atmosphere into homogeneous layers with constant properties. The “purely analytical” theories I to IV apply to a single isothermal atmospheric layer. What is proposed in this conclusion is actually a “hybrid analytical-numerical” approach that applies the analytical methods to each isothermal atmospheric layer, so as to obtain (iii) higher accuracy with a smaller number of steps over the altitude range of interest.
Considering the four analytical methods of calculation of rocket trajectories, each extends the typical textbook calculation of the effects of variable mass due to propellant consumption, by taking account of the dependence of the aerodynamic forces on the velocity and altitude. The methods I to IV significantly reduce the current large gap between the limited domain of validity of analytical solutions and the broader range of conditions covered by numerical methods. The issue here is not either numerical or analytical methods but rather the symbiotic complementarity of both: (i) the analytical methods make more explicit the effects of rocket, aerodynamic, and atmospheric parameters on the trajectory; and (ii) they assist in the interpretation of the results of numerical methods that can provide results in more complex situations, but whose interpretation may be less transparent. Thus, reducing the gap between analytical and numerical methods is beneficial for both approaches to the calculation of rocket trajectories.
The four methods I to IV of calculation of rocket trajectories were presented in the simplest one-dimensional case of vertical flight to show more clearly the effects of atmospheric conditions and rocket parameters, that appear in six quantities: (i) the dimensionless weight Equation (36a), (ii) thrust Equation (36b), and (iii) drag Equation (36c) parameters; (iv) the reference mass Equation (38a) that is the mass of a volume of air with the cross-sectional area of the rocket, the length of the atmospheric scale height at initial lift-off mass density, and can be made dimensionless dividing by the lift-off mass; (v) the reference time Equation (48b) that is the ratio of the lift-off mass to the rate of propellant consumption, and thus would be the burn time for zero residual mass, if the rocket mass was all propellant, and thus always exceeds the real burn time; (vi) the reference velocity Equation (48a) that corresponds to covering an atmospheric scale height in a reference time. Two distinct rockets with equal values of these six parameters would have the same trajectory specified by the methods I to IV. These methods can be extended from the one-dimensional case of a vertical climb to the practically more relevant two-dimensional case of a gravity turn and may be the subject of future work.

Author Contributions

Conceptualization, L.M.B.C.C.; Methodology, L.M.B.C.C.; Formal Analysis, L.M.B.C.C.; Resources, L.M.B.C.C. and P.J.S.G.; Writing—Original Draft Preparation, L.M.B.C.C. and P.J.S.G.; Writing—Review and Editing, L.M.B.C.C. and P.J.S.G.; Visualization, L.M.B.C.C. and P.J.S.G.; Project Administration, L.M.B.C.C. and P.J.S.G.

Funding

The work of the authors was supported by Fundação para a Ciência e a Tecnologia (FCT), through Institute of Mechanical Engineering (IDMEC), under the Associated Laboratory of Energy, Transport and Aeronautics (LAETA), project UID/EMS/50022/2013.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Culler, G.J.; Fried, B.D. Universal Gravity Turn Trajectories. J. Appl. Phys. 1957, 28, 672–676. [Google Scholar] [CrossRef]
  2. Thomson, W.T. Introduction to Space Dynamics, 2nd ed.; Dover Publications: New York, NY, USA, 1986. [Google Scholar]
  3. Miele, A. Flight Mechanics; Addison-Wesley: Boston, MA, USA, 1962; Volume 2. [Google Scholar]
  4. Connor, M.A. Gravity turn trajectories through the atmosphere. J. Spacecr. Rocket. 1966, 3, 1308–1311. [Google Scholar] [CrossRef]
  5. Sotto, E.D.; Teofilatto, P. Semi-Analytical Formulas for Launcher Performance Evaluation. J. Guid. Control Dyn. 2002, 25, 538–545. [Google Scholar] [CrossRef]
  6. Rutherford, D.E. Classical Mechanics, 2nd ed.; Oliver and Boyd: Edinburgh, Scotland, UK, 1967. [Google Scholar]
  7. Miele, A. Method of particular solutions for linear, two-point boundary-value problems. J. Optim. Theory Appl. 1968, 2, 260–273. [Google Scholar] [CrossRef]
  8. Miele, A.; Pritchard, R.E.; Damoulakis, J.N. Sequential gradient-restoration algorithm for optimal control problems. J. Optim. Theory Appl. 1970, 5, 235–282. [Google Scholar] [CrossRef]
  9. Von Mises, R. Theory of Flight; Dover: New York, NY, USA, 1959. [Google Scholar]
  10. Perkins, C.D.; Hage, R.E. Airplane Performance, Stability and Control; Wiley: Hoboken, NJ, USA, 1949. [Google Scholar]
  11. Milne-Thomson, L.M. Theoretical Aerodynamics; Dover: New York, NY, USA, 1958. [Google Scholar]
  12. Campos, L.M.B.C. Complex Analysis with Applications to Flows and Fields; CRC Press: Boca Raton, FL, USA, 2010. [Google Scholar]
  13. Campos, L.M.B.C. Mecânica Aplicada; Escolar Editora: Lisboa, Portugal, 2003. [Google Scholar]
  14. Ariane 5. Available online: https://en.wikipedia.org/wiki/Ariane_5 (accessed on 10 August 2018).
  15. Ariane 5—VA226—Launch Profile. Available online: https://spaceflight101.com/ariane-5-va226/ariane-5-va226-launch-profile/ (accessed on 10 August 2018).
Figure 1. Forces on a rocket in flight. See text for the definition of the symbols.
Figure 1. Forces on a rocket in flight. See text for the definition of the symbols.
Aerospace 05 00088 g001
Figure 2. Approximation to the trajectory of a rocket by homogeneous altitude layers.
Figure 2. Approximation to the trajectory of a rocket by homogeneous altitude layers.
Aerospace 05 00088 g002
Figure 3. Approximation to the trajectory of a rocket by arcs of constant curvature.
Figure 3. Approximation to the trajectory of a rocket by arcs of constant curvature.
Aerospace 05 00088 g003
Figure 4. Approximation to the trajectory of a rocket by tangents with a constant flight path angle.
Figure 4. Approximation to the trajectory of a rocket by tangents with a constant flight path angle.
Aerospace 05 00088 g004
Table 1. Effects included in the calculation of rocket trajectories in the literature [1,2,3,4,5] and in this work.
Table 1. Effects included in the calculation of rocket trajectories in the literature [1,2,3,4,5] and in this work.
NoEffectPresent Work[1][2][3][4][5]
1rocket mass as function of time
2drag included
3lift included
4aerodynamic forces depend on velocity
5aerodynamic forces depend on altitude
6flight path angle variable
7non-zero angle-of-attack
8thrust vector angle to flight path
9multistage
10wind
Table 2. Input data for calculations of rocket trajectories by methods I in Section 3, method II in Section 4 and method III in Section 5.1 and Section 5.2.
Table 2. Input data for calculations of rocket trajectories by methods I in Section 3, method II in Section 4 and method III in Section 5.1 and Section 5.2.
SymbolMeaningValueUnit
Environment
gacceleration of gravity9.81 m   s 1
ρ 0 sea level mass density1.225 k g   m 3
atmospheric scale height 2.6 × 10 4 m
Propulsion
Tthrust 1.555 × 10 7 k g   m   s 2
cpropellant flow rate 2.029 × 10 3 k g   s 1
m 0 m 1 propellant mass 2.841 × 10 5 k g
t 1 t 0 burn time 1.4 × 10 2 s
Aerodynamics
m 0 initial mass 7.77 × 10 5 k g
Scross-sectional area3.76 × 10 m 2
C D drag coefficient0.15*
Calculated parameters
m * reference mass 1.198 × 10 6 k g
aweight parameter5.533 × 10*
bthrust parameter 1.129 × 10 2 *
fdrag parameter 1.156 × 10 1 *
t * reference time 3.829 × 10 2 s
v * reference velocity6.790 × 10 m   s 1
* Dimensionless parameter.
Table 3. Calculation of powered trajectories with four discretizations for method I in Section 6.1 and method III in Section 6.2.
Table 3. Calculation of powered trajectories with four discretizations for method I in Section 6.1 and method III in Section 6.2.
MethodIIII
Coefficients of series
first A 0 = 1.00 D 0 = 0
second A 1 = 0.00 D 1 = 0
third A 2 = 28.77 D 2 = 5.10 m   s 2
fourth A 3 = 18.81 D 3 = 8.71 × 10 3 m   s 3
fifth A 4 = 391.4 N/A
Altitude versus time ( m )
t = 35   s 0.646 × 10 4 0.662 × 10 4
t = 70   s 2.396 × 10 4 2.798 × 10 4
t = 105   s 4.554 × 10 4 6.633 × 10 4
t = 140   s 6.620 × 10 4 1.239 × 10 5
Velocity versus time ( m   s 1 )
t = 35   s 3.668 × 10 2 3.891 × 10 2
t = 70   s 5.935 × 10 2 8.422 × 10 2
t = 105   s 6.161 × 10 2 1.359 × 10 3
t = 140   s 5.593 × 10 2 1.941 × 10 3
Table 4. Calculation of ballistic trajectory by method IV in Section 5.3, using initial data from methods I or II in Section 6.3.
Table 4. Calculation of ballistic trajectory by method IV in Section 5.3, using initial data from methods I or II in Section 6.3.
MethodI + IVIII + IVUnits
Burnout data ( m 1 = 4.929 × 10 5   k g )
z 1 6.610 × 10 4 1.239 × 10 5 m
z ˙ 1 5.593 × 10 2 1.941 × 10 3 m   s 1
ρ 1 = ρ 0 e z / 9.602 × 10 2 1.044 × 10 2 k g   m 3
Calculated values
m * * = ρ 1 S 9.390 × 10 4 1.021 × 10 4 k g
ϑ = C D m * * / m 1 2.857 × 10 2 3.107 × 10 3 *
E ¯ 0 = z ˙ 1 2 / ( 2 g ) 6.132 × 10 1 7.382*
E 0 = E ¯ 0 e ϑ 5.959 × 10 1 7.359*
Apogee of trajectory ( m 1 = 4.929 × 10 5 k g )
X = E 0 / ( 1 ϑ ) 6.134 × 10 1 7.382*
z 2 z 1 = X 1.595 × 10 4 1.919 × 10 5 m
z 2 8.215 × 10 4 3.158 × 10 5 m
* Dimensionless quantities.

Share and Cite

MDPI and ACS Style

Campos, L.M.B.C.; Gil, P.J.S. On Four New Methods of Analytical Calculation of Rocket Trajectories. Aerospace 2018, 5, 88. https://doi.org/10.3390/aerospace5030088

AMA Style

Campos LMBC, Gil PJS. On Four New Methods of Analytical Calculation of Rocket Trajectories. Aerospace. 2018; 5(3):88. https://doi.org/10.3390/aerospace5030088

Chicago/Turabian Style

Campos, Luís M. B. C., and Paulo J. S. Gil. 2018. "On Four New Methods of Analytical Calculation of Rocket Trajectories" Aerospace 5, no. 3: 88. https://doi.org/10.3390/aerospace5030088

APA Style

Campos, L. M. B. C., & Gil, P. J. S. (2018). On Four New Methods of Analytical Calculation of Rocket Trajectories. Aerospace, 5(3), 88. https://doi.org/10.3390/aerospace5030088

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop