Trajectory Optimization of Launch Vehicles Using Object-oriented Programming

ABSTRACT:

The aim of this study is to model launch vehicles with focus on 3-DOF trajectory optimization using a modular approach. Despite the large number of operational launch vehicles, they usually consist of basic components and subsystems. In other words, a launch vehicle is an assembly of stages, which in turn is divided into propellant system and engine, and the engine is an assembly of basic components such as pumps, turbines, combustion chamber, and nozzle. To allow future extension and reuse of the codes, a modular structure using object-oriented programming is used. Two formulations of state equations of the trajectory and two optimization methods are described. The launch vehicle performance will be measured by payload mass for a given mission. The simulations of the VLS-1, Ariane 5 and VLS-Alfa were performed and showed good agreement with the literature.

KEYWORDS:
Launch vehicle; Trajectory; Optimization; Object-oriented programming

INTRODUCTION

Due to the inherent complexity of a launch vehicle, its design is traditionally divided into multiple disciplines, such as trajectory optimization, propulsion, aerodynamics, and mass budget. Despite the large number of operational launch vehicles, they usually consist of basic components. In other words, a launch vehicle is an assembly of stages, which in turn is divided into propellant system and engine, and the engine is an assembly of basic components such as pumps, turbines, combustion chamber, and nozzle (for a liquid rocket stage). Space launch systems are composed by a large number of components grouped into a hierarchy of subsystems. The performance of the vehicle depends on the individual performance of each of the subsystems, which in turn depend on material properties and design parameters. Changes in design parameters are propagated throughout the cluster hierarchy of subsystems and components, flight trajectory and payload capability (Hinckel 1995).

In order to get the best performance of a given launch vehicle, and consequently, to make the access to space less costly, trajectory optimization techniques have been for decades a subject of intense research. Trajectory optimization can be categorized basically into direct and indirect methods. Betts (1998) and Rao (2009) have made a comprehensive discussion about both methods. The reason for this method to be called "indirect" comes from the strategy to convert the original optimal control problem into a boundary-value problem. The most common indirect methods found in the literature are the shooting method, the multiple-shooting method, and collocation methods (Brown et al. 1969; Teren and Spurlock 1966; Miele and Wang 2003). Presumably because of the possibility of solving complex problems with a minimum effort of mathematical analysis, the direct method is the one chosen by most researchers (Hargraves and Paris 1987; Seywald 1994; Herman and Conway 1996; Balesdent 2011). One of the most popular software used extensively in many publications is called POST (Program to Optimize Simulated Trajectories) (Brauer et al. 1977). POST is a commercial program that has been used successfully to solve a wide variety of atmospheric ascent and reentry problems, as well as exoatmospheric orbital transfer problems. In the framework of this method, the problem is characterized by a set of parameters that define the control law. This problem is a typical Non Linear Programming Problem (NLP) and can be solved using classical Gradient-based methods (deterministic methods) such as Sequential Quadratic Program (SQP) or by heuristic methods. According to Betts (1998), heuristic optimization algorithms are not computationally competitive with gradient methods. Even though presumably due to ease of implementation without a detailed understanding of the system, in the last two decades lots of papers using Particle Swarm Optimization (PSO), genetic algorithms (GA), among others were applied to solve trajectory optimization problems. As for indirect methods, the direct methods can be categorized in direct (multiple) shooting or collocation. In the case where only the control variables are adjusted by a function, the method is called a shooting method. When both the state and control are parameterized, the method is called a collocation method. A well-known software developed by the University of Stuttgart, which addresses the direct collocation method, is the AeroSpace Trajectory Optimization Software (ASTOS).

For either direct or indirect approaches, perhaps the most important benefit gained from a multiple shooting formulation compared to its precursor (single shooting) is enhanced robustness. To take advantage of both methods previously described, a hybrid method can also be considered (von Stryk and Bulirsch 1992; Pontani and Teofilatto 2014; Gath and Calise 2001). The idea behind this approach is to divide the flight trajectory into two distinct phases, namely atmospheric and exoatmospheric phases, applying the direct method in the first phase and indirect method in the second one. Here, exoatmospheric phase means that the vehicle is virtually in vacuum space, i.e., the aerodynamic effects can be ignored. Pontani and Teofilatto (2014) proposed a simple method to evaluate the performance of multistage launch vehicles for given structural data, aerodynamic and propulsive parameters.

The purpose of this work is to develop a tool, which can be easily reused and extended, to model and simulate a launch vehicle with focus on trajectory. Thus, the tool is intended to be capable of modeling and optimizing flight trajectory until orbit injection. Two mathematical models for determination of performance of a launch vehicle will be described and discussed. The launch vehicle performance will be measured by payload mass for a given mission.

TRAJECTORY MODELING

Since the dawn of the space era launch vehicles are responsible to put satellites into orbit. This makes the cost of a satellite strongly related to the performance of the launch vehicle, which in turn depends on the trajectory profile.

ATMOSPHERE MODEL

The atmosphere can be seen as a layer of gases attached to the surface of the Earth by gravitational attraction. This work makes use of the standard atmosphere, which is modeled as adjacent layers of gases in which temperature depends on the altitude. In different layers, the temperature can be modeled as linear function of the altitude.

AERODYNAMICS

During the flight, a launch vehicle needs to cross the atmosphere in which reacts to the vehicle motion by means of aerodynamics forces. The drag force arises due to friction between the body and the fluid (Eq. 1):

(1) $\mathit{D}=\frac{1}{2}\mathit{\rho }\left(\mathit{r}\right){\mathit{SC}}_{\mathit{D}}{\mathit{V}}^2$

where: D = drag force (N); ρ = density of the working fluid (kg/m³); S = reference area (m²); C_D = drag coefficient (-); and V = air speed (m/s).

The lift is a reaction force to the angle of attack (Eq. 2):

(2) $\mathit{L}=\frac{1}{2}\mathit{\rho }\left(\mathit{r}\right){\mathit{SC}}_{\mathit{L}}{\mathit{V}}^2$

where: L = lift force (N); S = reference area (m²); C_L = lift coefficient (-); and V = air speed (m/s).

There are many methods to estimate the aerodynamic coefficients in the literature, i.e., the well-known tool Missile DATCOM, interpolation of available data from a given vehicle, closed formulas that consider contributions of shock wave in the rocket nose, body friction and base pressure, or even constant value for certain phases of the flight. Except for interpolation from real flight data, the methods to estimate the drag coefficient are quite inaccurate. Fortunately, because the essential acceleration phase begins in the exoatmospheric phase, usually the drag has little influence on launcher performance (Schlingloff 2005).

EQUATIONS OF THE TRANSLATIONAL MOTION

The modeling of the trajectory of a launch vehicle is usually performed by means of two reference frames (one with origin on the Earth center and the other one moving with the vehicle) and considerations or idealizations according to the requirements of the mission. To model the translational motion, the vehicle can be treated as a particle, ignoring the size and mass distribution. In modeling the rotational motion, the vehicle can be considered a rigid body, reducing the degrees of freedom from infinity (flexible body case) to just six (Tewari 2007). Strictly speaking, a launch vehicle is far from being considered a rigid body. Mass is continuously expelled due to combustion of the propellants and for liquid rocket engines there is still the sloshing, which is the movement of fluid within the tanks and pipes and rotating equipment such as turbines and pumps. Especially for large launch vehicles, the deflection of the structure should be considered as well (Cornelisse et al. 1979). However, this research is focused on the reference trajectory, thus treatment of the translational motion is sufficient to fulfill this task. Two mathematical models for the trajectory are presented in the following sections. The first one was taken from Schlingloff (2005) and the second one from Tewari (2007).

First Formulation

In this formulation, the spherical celestial (inertial) coordinates and a moving coordinates in the orbit plane were considered. Both reference frames have the origin on the Earth center (Fig. 1). The vector of state variables is conveniently chosen as y(t) = [r(t) u(t) v(t) Ω(t) ι(t) ω(t)]^T, where u = vertical velocity (m/s); v = horizontal velocity (m/s); where: Ω = right ascension of the ascending node (rad) and ω = argument of periapsis (rad). Thus, the system of equation can be given as (Eqs. 3-9):

Figure 1
Reference frames.

(3) $\dot{\mathit{r}}=\mathit{u}$

(4) $\dot{\mathit{u}}=\frac{{\mathit{v}}^2}{\mathit{r}}-\frac{\mathit{\mu }}{{\mathit{r}}^2}+\frac{\mathit{F}}{\mathit{m}}\sin \mathit{\beta }\cos \mathit{\delta }+\frac{{\mathit{D}}_{\mathit{x}}}{\mathit{m}}+\frac{{\mathit{L}}_{\mathit{x}}}{\mathit{m}}$

(5) $\dot{\mathit{v}}=-\frac{\mathit{\mu }\mathit{v}}{\mathit{r}}+\frac{\mathit{F}}{\mathit{m}}\cos \mathit{\beta }\cos \mathit{\delta }+\frac{{\mathit{D}}_{\mathit{y}}}{\mathit{m}}+\frac{{\mathit{L}}_{\mathit{y}}}{\mathit{m}}$

(6) $\dot{\Omega }={\mathit{\omega }}_{\mathit{x}}\frac{\sin \mathit{\omega }}{\sin \mathit{\iota }}$

(7) $\mathit{i}={\mathit{\omega }}_{\mathit{x}}\cos \mathit{\omega }$

(8) $\dot{\mathit{\omega }}=\frac{\mathit{v}}{\mathit{r}}-{\mathit{\omega }}_{\mathit{x}}\sin \mathit{\omega }\cot \mathit{\iota }$

(9) ${\mathit{\omega }}_{\mathit{x}}=\frac{\mathit{F}\sin \mathit{\delta }+{\mathit{D}}_{\mathit{z}}+{\mathit{L}}_{\mathit{z}}}{\mathit{v}}$

where: β = thrust angle in flight plane (rad); δ = thrust angle out of flight plane (rad); and F = thrust force (N); ω_x = inclination change (rad).

Equations 3 to 5 of the system of differential equations are the dynamic equations of motion and Eqs. 6 to 9 are the kinematic equations. The dynamic equations are derived by application of the Newton's Second Law resolved into components of the moving system. The kinematic equations are deducted into two steps: representation of the rotation velocity of the vehicle in a vector form, and applying the Euler angles (Schlingloff 2005).

Second Formulation

The reference frames adopted in this modeling are the planet-fixed reference (SXYZ) frame and the local horizontal frame (oxyz), both are non-inertial (Fig. 2).

Figure 2
Planet-fixed and local horizon frames for atmospheric flight (adapted from Tewari 2007).

From Fig. 2, the relative velocity v and the local velocity of the local horizontal frame (oxyz) relative to the planet-centered rotating frame (SXYZ) can be expressed as (Eqs. 10 and 11):

(10) $\mathbf{v}\left(\mathit{v},\mathit{y},\mathit{\zeta }\right)=\mathit{v}\left(\sin \mathit{\gamma }\mathbf{i}+\cos \mathit{\gamma }\sin \mathit{\zeta }\mathbf{j}+\cos \mathit{\gamma }\cos \mathit{\zeta }\mathbf{k}\right)$

(11) $\Omega =\dot{\mathit{\xi }}\mathbf{K}-\dot{\mathit{\varphi }}\mathbf{j}$

where: γ = flight path angle (rad); ζ = heading angle (= π/2 - A); and ξ = longitude (rad); and A = azimuth (rad).

With a convenient rotation matrix, Eq. 11 can be written only in terms of axes of the body as (Eq. 12):

(12) $\mathbf{\Omega }=\dot{\mathit{\xi }}\sin \mathit{\varphi }\mathbf{i}-\dot{\mathit{\varphi }}\mathbf{j}+\dot{\mathit{\xi }}\cos \mathit{\varphi }\mathbf{k}$

The relative velocity can also be expressed as (Eqs. 13 and 14):

(13) $\mathbf{v}=\dot{\mathit{r}}\mathbf{i}+\mathbf{\Omega }\times \left(\mathit{r}\mathbf{i}\right)$

(14) $\mathbf{v}=\dot{\mathit{r}}\mathbf{i}+\mathit{r}\dot{\mathit{\xi }}\cos \mathit{\varphi }\mathbf{j}+\mathit{r}\dot{\mathit{\varphi }}\mathbf{k}$

Comparing Eqs. 13 and 14 we finally obtain the kinematic equations of motion (Eqs. 15-17):

(15) $\dot{\mathit{r}}=\mathit{v}\sin \mathit{\gamma }$

(16) $\dot{\mathit{\xi }}=\frac{\mathit{v}\cos \mathit{\gamma }\cos \mathit{\zeta }}{\mathit{r}\cos \mathit{\varphi }}$

(17) $\dot{\mathit{\varphi }}=\frac{\mathit{v}\cos \mathit{\gamma }\cos \mathit{\zeta }}{\mathit{r}}$

To derive the dynamic equations Newton's Second Law must be introduced (Eq. 18):

(18) $\mathbf{f}=\mathit{m}{\mathbf{a}}_{\mathbf{I}}=\mathit{m}\frac{\mathit{d}{\mathbf{v}}_{\mathbf{I}}}{\mathit{dt}}$

Choosing the wind axes to express the forces on the body and doing the appropriate transformation to perform a_I in the wind axes, the remaining equations to model the translational motion are obtained (Eqs. 19-21):

(19) $\dot{\mathit{\gamma }}=\frac{\mathit{F}\sin {\mathit{\alpha }}_{\mathit{T}}}{\mathit{mv}}+\left(\frac{\mathit{v}}{\mathit{r}}-\frac{{\mathit{\mu }}_{\mathit{E}}}{{\mathit{r}}^2\mathit{v}}\right)\cos \mathit{\gamma }+\frac{\mathit{L}}{\mathit{mv}}+\cos \mathit{\varphi }\left[2{\mathit{\omega }}_{\mathit{E}}\cos \mathit{\zeta }+\frac{{\mathit{\omega }}_{\mathit{E}}^2\mathit{r}}{\mathit{v}}\left(\cos \mathit{\varphi }\cos \mathit{\gamma }+\sin \mathit{\varphi }\sin \mathit{\gamma }\sin \mathit{\zeta }\right)\right]$

(20) $\dot{\mathit{v}}=\frac{\mathit{F}\cos {\mathit{\alpha }}_{\mathit{T}}}{\mathit{m}}-\frac{{\mathit{\mu }}_{\mathit{E}}}{{\mathit{r}}^2}\sin \mathit{\gamma }-\frac{\mathit{D}}{\mathit{m}}+{\mathit{\omega }}_{\mathit{E}}^2\mathit{r}\cos \mathit{\varphi }\left(\cos \mathit{\varphi }\sin \mathit{\gamma }-\sin \mathit{\varphi }\cos \mathit{\gamma }\sin \mathit{\zeta }\right)$

(21) $\dot{\mathit{\zeta }}=-\frac{\mathit{v}}{\mathit{r}}\tan \mathit{\varphi }\cos \mathit{\gamma }\cos \mathit{\zeta }+2{\mathit{\omega }}_{\mathit{E}}\cos \mathit{\varphi }\tan \mathit{\gamma }\sin \mathit{\zeta }-\frac{{\mathit{\omega }}_{\mathit{E}}^2\mathit{r}}{\mathit{v}\cos \mathit{\gamma }}\sin \mathit{\varphi }\cos \mathit{\varphi }\cos \mathit{\zeta }-2{\mathit{\omega }}_{\mathit{E}}\sin \mathit{\varphi }$

where: ω_E = Earth rotation (rad/s).

Equations 15-17 are the kinematical equations of motion and Eqs. 19-21 are the dynamic equations. With the integration of the system of differential equations, the vector position and the vector velocity of the vehicle can be determined by the following equations (Eqs. 22 and 23):

(22) $\mathbf{r}\left(\mathit{r},\mathit{\varphi },\mathit{\xi }\right)=\mathit{r}\left(\cos \mathit{\varphi }\cos \mathit{\xi }\mathbf{I}+\cos \mathit{\varphi }\sin \mathit{\xi }\mathbf{J}+\sin \mathit{\varphi }\mathbf{K}\right)$

(23) $\mathbf{v}\left(\mathit{v},\mathit{\gamma },\mathit{\zeta }\right)=\mathit{v}\left(\sin \mathit{\gamma }\mathbf{i}+\cos \mathit{\gamma }\sin \mathit{\zeta }\mathbf{j}+\cos \mathit{\gamma }\cos \mathit{\zeta }\mathbf{k}\right)$

It's known that if one has an inertial vector position and a velocity vector of a given body in orbit, the orbital elements (or Keplerian elements) can be readily determined. Thus, to get the orbital elements, it is necessary to perform an appropriate matrix rotation to obtain the desired inertial vectors.

OPTIMIZATION

In order to obtain the maximal payload capacity of a given launch vehicle, and consequently, make the access to space less expensive, trajectory optimization techniques have been for decades a subject of intense research. The trajectory optimization can be categorized basically into direct and indirect methods. To take advantage of both methods, a combination of both techniques can also be done, i.e., a hybrid method can also be considered.

METHODOLOGY

In this section the techniques used to solve the trajectory optimization problem are presented. The first approach was based on Silva (1995). The second one describes a hybrid algorithm that merges the direct and indirect methods.

FIRST APPROACH - DIRECT METHOD

The method applied within the framework of this approach is based on Silva (1995). A polynomial control function is used to model the flight profile. Four parameters, namely the coast time duration t_coast (if it is applied) and three parameters of the polynomial control function, are optimized in order to get the maximum payload mass. A code from Jacob (1972) written in FORTRAN is transcript to C++ language and adapted to solve the problem (Eq. 24).

(24) $\begin{array}{c}\mathit{\beta }=\left\{\begin{array}{c}\mathit{\pi }/180,\mathrm{if}\mathit{t}\le {\mathit{t}}_{\mathit{v}}\\ {\mathit{b}}_0-{\mathit{b}}_1\left(\mathit{t}-{\mathit{t}}_{\mathit{v}}\right)+{\mathit{b}}_2{\left(\mathit{t}-{\mathit{t}}_{\mathit{v}}\right)}^2,\mathrm{if}{\mathit{t}}_{\mathit{v}}<\mathit{t}\le {\mathit{t}}_{\mathit{b}1}\\ {\mathit{b}}_3-{\mathit{b}}_4\left(\mathit{t}-{\mathit{t}}_{\mathit{f}1}\right)+{\mathit{b}}_5{\left(\mathit{t}-{\mathit{t}}_{\mathit{b}1}\right)}^2,\mathrm{if}{\mathit{t}}_{\mathit{b}1}<\mathit{t}\le {\mathit{t}}_{\mathit{bf}}\end{array}\right.\hfill \\ \mathrm{where}:{\mathit{b}}_0=\frac{\mathit{\pi }}{2};{\mathit{b}}_1=\frac{{\mathit{\beta }}_1\mathit{\pi }}{180\left({\mathit{t}}_{\mathit{b}1}-{\mathit{t}}_{\mathit{v}}\right)};{\mathit{b}}_2=-\frac{{\mathit{b}}_0-{\mathit{\beta }}_2\mathit{\pi }/180-{\mathit{b}}_1\left({\mathit{t}}_{\mathit{b}1}-{\mathit{t}}_{\mathit{v}}\right)}{{\left({\mathit{t}}_{\mathit{b}1}-{\mathit{t}}_{\mathit{v}}\right)}^2};{\mathit{b}}_3=\frac{{\mathit{\beta }}_2\mathit{\pi }}{180};\hfill \\ {\mathit{b}}_4=\frac{{\mathit{\beta }}_3\mathit{\pi }}{180\left({\mathit{t}}_{\mathit{bf}}-{\mathit{t}}_{\mathit{b}1}\right)};{\mathit{b}}_5=\frac{{\mathit{b}}_3-{\mathit{b}}_4\left({\mathit{t}}_{\mathit{bf}}-{\mathit{t}}_{\mathit{b}1}\right)}{\left({\mathit{t}}_{\mathit{bf}}-{\mathit{t}}_{\mathit{b}1}\right)};\hfill \end{array}$

and β₁, β₂, β₃ = set of optimization control parameters (-).

The trajectory optimization problem can be formulated as (Eq. 25):

(25) $\mathrm{Find}\mathbf{X}=\left[\begin{array}{c}{\mathit{t}}_{\mathit{coast}}\\ {\mathit{\beta }}_1\\ {\mathit{\beta }}_2\\ {\mathit{\beta }}_3\end{array}\right]$

which maximize the pay load mass subjects to the constraints at orbit injection. Thus, for the first formulation of equations of motion (Eqs. 26-28):

(26) $\mathit{r}={\mathit{R}}_{\mathit{e}}+{\mathit{h}}_{\mathit{f}}$

(27) $\mathit{u}=0$

(28) $\mathit{v}=\sqrt{\frac{\mathit{\mu }}{\mathit{r}}}$

where: R_e = equatorial radius of the Earth (km).

For the second approach (Eqs. 29-31):

(29) $\mathit{r}={\mathit{R}}_{\mathit{e}}+{\mathit{h}}_{\mathit{f}}$

(30) $\mathit{\gamma }=0$

(31) ${\mathit{v}}_{\mathit{I}}=\sqrt{{\mathit{v}}^2+{\left({\mathit{\omega }}_{\mathit{E}}\mathit{r}\cos \mathit{\varphi }\right)}^2+2\mathit{v}{\mathit{\omega }}_{\mathit{E}}\mathit{r}\cos \mathit{\varphi }\cos \mathit{\gamma }\cos \mathit{\zeta }}=\sqrt{\frac{\mathit{\mu }}{\mathit{r}}}$

SECOND APPROACH - HYBRID METHOD

The idea behind the strategy is to divide the flight trajectory into two distinct phases, one while the vehicle ascents the dense atmosphere and the other one when the vehicle is virtually in vacuum space, i.e., where aerodynamic effects can be ignored. In the aforementioned sections it was stated that the sensibility of the indirect shooting or even multiple shooting method depends on the initial guess, however, it is possible to compute the initial co-states variable for optimal thrust arcs in vacuum fairly easy using almost arbitrary initial guess when direct and indirect methods are merged.

STAGES BEFORE COASTING - DIRECT METHOD

To fulfill this task, the method used in the first method is applied here. In other words, the same control variables will be obtained until the beginning of the coast arc. Furthermore, this step gives the gross lift-off mass (GLOW) of the vehicle.

UPPER STAGE TRAJECTORY - INDIRECT METHOD

In this phase, the vehicle is assumed to be out of the denser layers of the atmosphere, so that the aerodynamics forces can be neglected. Thus, the trajectory of the upper stage is accomplished in the orbit plane. The upper stage flight is divided into two phases: a coast arc and a thrust arc. The equations of motion are presented below (Eqs. 32-34) (Schlingloff 1987).

(32) $\dot{\mathit{r}}=\mathit{u}$

(33) $\dot{\mathit{u}}=\frac{{\mathit{v}}^2}{\mathit{r}}-\frac{\mathit{\mu }}{{\mathit{r}}^2}+\frac{\mathit{F}}{\mathit{m}}\sin \mathit{\beta }$

(34) $\dot{\mathit{v}}=-\frac{\mathit{\mu }\mathit{v}}{\mathit{r}}+\frac{\mathit{F}}{\mathit{m}}\cos \mathit{\beta }$

Applying the Lagrange method, the Hamiltonian H is constructed as (Eq. 35):

(35) $\mathit{H}={\mathit{\lambda }}_{\mathit{u}}\left(\frac{{\mathit{v}}^2}{\mathit{r}}-\frac{\mathit{\mu }}{{\mathit{r}}^2}+\frac{\mathit{F}}{\mathit{m}}\sin \mathit{\beta }\right)+{\mathit{\lambda }}_{\mathit{v}}\left(-\frac{\mathit{\mu }\mathit{v}}{\mathit{r}}+\frac{\mathit{F}}{\mathit{m}}\cos \mathit{\beta }\right)+{\mathit{\lambda }}_{\mathit{r}}\mathit{u}$

The adjoint equations for the co-state (Euler-Lagrange equations) are (Eqs. 36-38):

(36) ${\dot{\mathit{\lambda }}}_{\mathit{u}}=-\frac{\partial \mathit{H}}{\partial \mathit{u}}=\frac{\mathit{v}}{\mathit{r}}{\mathit{\lambda }}_{\mathit{v}}-{\mathit{\lambda }}_{\mathit{r}}$

(37) ${\dot{\mathit{\lambda }}}_{\mathit{v}}=-\frac{\partial \mathit{H}}{\partial \mathit{v}}=-2\frac{\mathit{v}}{\mathit{r}}{\mathit{\lambda }}_{\mathit{u}}+\frac{\mathit{u}}{\mathit{r}}{\mathit{\lambda }}_{\mathit{v}}$

(38) ${\dot{\mathit{\lambda }}}_{\mathit{u}}=-\frac{\partial \mathit{H}}{\partial \mathit{u}}=\left(\frac{{\mathit{v}}^2}{{\mathit{r}}^2}-2\frac{\mathit{\mu }}{{\mathit{r}}^3}\right){\mathit{\lambda }}_{\mathit{u}}-\frac{\mathit{uv}}{{\mathit{r}}^2}{\mathit{\lambda }}_{\mathit{v}}$

Using the Pontryagin minimum principle, the optimal thrust angle control can be expressed as a function of the co-states (Eq. 39):

(39) $\frac{\partial \mathit{H}}{\partial \mathit{\beta }}=\frac{\mathit{F}}{\mathit{m}}\left({\mathit{\lambda }}_{\mathit{u}}\cos \mathit{\beta }-{\mathit{\lambda }}_{\mathit{v}}\sin \mathit{\beta }\right)=0$

and finally (Eq. 40):

(40) $\tan \mathit{\beta }=\frac{{\mathit{\lambda }}_{\mathit{u}}}{{\mathit{\lambda }}_{\mathit{v}}}$

Thus, to get the optimal trajectory, Eqs. 32-34 along with Eqs. 36-38 should be integrated. However, as the Lagrange multipliers have no physical meaning and the flight path depends very sensitively on the initial guesses, it is difficult to solve these equations. Schlingloff (1987) developed an analytical method to eliminate the Lagrange multipliers getting formulas that can be represented by a smaller number of variables. Thus, defining z (Eq. 41):

(41) $\mathit{z}=\frac{{\mathit{r}}^2\dot{\mathit{\beta }}-\mathit{rv}}{{\cos }^2\mathit{\beta }}$

where: z = variable with no physical meaning.

Schlingloff (1987) got alternative differential equations to represent the control law (Eqs. 42 and 43):

(42) $\tan \dot{\mathit{\beta }}=\frac{\mathit{z}}{{\mathit{r}}^2}+\frac{\mathit{v}}{\mathit{r}}\left(\tan {\mathit{\beta }}^2+1\right)$

(43) $\dot{\mathit{z}}=\frac{\tan \mathit{\beta }}{\mathit{r}}\left(4\mathit{vz}+3\mathit{y}\right)$

Thus, to get the new system of first-order differential equations, the control Eqs. 42 and 43 must be joined to the equations of motion Eqs. 32-34. To integrate this system, since the initial condition for the state equations are fixed, just the initial guesses to the control equations must be set. To take into account the coast, the optimization problem can be stated as:

Find X = [β₀ z₀ t_coast]^T which minimize the propellant mass of the upper stage F(X) = m_prop. It implies to maximize the payload mass that can be injected into the desired orbit. The constraints at orbit injection are given by Eqs. 26-28. To solve this problem we can use heuristic methods such as particle swarm optimization (PSO).

PROGRAM STRUCTURE

As aforementioned stated, a modular approach using object-oriented programming (OOP) is chosen and to allow a better visualization of the codes it is used UML diagrams. UML diagrams are used to visualize the code and the communication between objects enabling a high degree of abstraction. Figure 3 presents a UML diagram for a specific launch vehicle, namely the VLS-alfa launch vehicle. From the diagram we can see some parameters and functions of each component and the relationship between them. In order to make the diagram clearer, some parameters and functions are omitted.

Figure 3
Schematic of a UML diagram of the VLS-Alfa launch vehicle.

RESULTS

To verify the trajectory program, two launch vehicles are considered, namely: the Brazilian VLS-1 and the European Ariane 5. Both mathematical modeling of the ascent trajectory are used in order to verify the concordance between them. The trajectory optimization will be performed using direct method.

VLS-1 LAUNCH VEHICLE

The under development VLS-1 is the future Brazilian satellite launch vehicle. Its development started in 1984, however due to technical problems, the vehicle could not be qualified up to now. VLS-1 is composed of four solid stages. The first stage is equipped with four solid boosters S43 (Fig. 4). The vehicle is designed to perform a non-powered cost arc between third and upper stages. Key parameters of the vehicle used in the simulation are given in Fig. 4. The mission is to launch a satellite into a reference circular orbit of 500 km of altitude from the Alcântara Launch Center (2°22'39.52" S, 44°23'57.71" W).

Figure 4
VLS-1: key parameters.

Figures 5 and 6 show the velocity and altitude profiles with the powered phase described by the red curve, and the non-powered phase by the blue curve. A path constraint (dynamic pressure) of the flight is shown in Fig. 7. The trajectory profile of the VLS-1 was verified with the SKYNAV tool (Schlingloff 1991) and presented a good agreement. Tables 1 and 2 summarize the corner instants between flight phases and Table 3 presents the control parameters obtained by optimization subroutine. As expected, the maximum payload mass found for both formulations are very close.

Figure 5
Velocity profile of the VLS-1 launch vehicle using (a) first formulation of state equations and (b) the second formulation.

Figure 6
Altitude profile of the VLS-1 launch vehicle using (a) first and (b) the second formulation of state equations.

Figure 7
Dynamic pressure profile of the VLS-1 launch vehicle using (a) first and (b) the second formulation of state equations.

Table 1
Values for state variables at start, inter-stage and end instants.

t (s)	h (km)	u (m/s)	v (m/s)	Ω (deg)	ι (deg)	ω (deg)
0.0	0.0	0.0	421.5	90.0	-2.4	-44.4
62.8	23.5	763.5	1390.0	90.4	-2.4	-63.3
124.9	75.4	1023.9	2729.1	91.5	-2.4	-64.1
183.2	154.5	1927.9	4947.7	93.4	-2.4	-64.1
530.9	498.1	69.8	4700.4	107.5	-2.4	-64.1
605.4	499.8	1.2	7612.6	111.2	-2.4	-64.1

Table 2
Values for state variables at start, inter-stage and end instants.

t (s)	h (km)	v (m/s)	ξ (deg)	φ (deg)	γ (deg)	ζ (deg)
0.0	0.0	0.0	-44.4	-2.4	90	0
62.8	23.4	1211.0	-44.2	-2.4	38.9	0.4
124.9	75.2	2493.5	-43.4	-2.4	24.2	0.4
183.2	154.4	4881.4	-41.8	-2.4	23.3	0.5
530.9	498.0	4213.5	-29.1	-2.2	1.0	1.1
605.5	499.9	7116.8	-25.7	-2.1	0.0	1.3

Table 3
Optimized control parameters.

#	mp (kg)	tcoast (s)	β1	β2	β3
1st formulation	270.5	347.71	72.57	31.47	-0.08
2nd formulation	269.40	347.48	72.86	31.61	-0.07

Comparing the results from the Tables 1 to 3 and Figs. 5 to 7, we can verify that indeed both mathematical modeling are equivalent. Although the method presented in this work gives sub-optimal trajectory, for the purpose of a preliminary analysis, this method is sufficiently accurate.

ARIANE 5 LAUNCH VEHICLE

Built under supervision of European Space Agency (ESA), Ariane 5 is a European launch vehicle that is part of the Ariane rocket family. The vehicle is used to deliver payload into low Earth orbit (LEO) and geostationary transfer orbit (GTO). Within the framework of this work, the mission is to launch a satellite from Kourou to a low Earth orbit (LEO) of 200 km of altitude. The key parameters of the vehicle are given in Table 4.

Table 4
Key parameters: Ariane 5 launch vehicle (Schlingloff 1987).

Ariane 5	mprop (ton)	ms (ton)	Isp (s)	tb (s)	CD*S (m2)
1st Stage	155	15	342/439	651.756	18.32
2nd Stage	450	80	245/279	146.84	30.396
3rd Stage	2.4	0.7	320	384	0.0

The trajectory can be divided into two main phases. To take the vehicle from the ground and minimize the gravitational losses, the first phase is powered by two solid boosters and by the core stage using the propellants LOX/LH2. After 146.84 seconds the booster stages are decoupled from the vehicle and the motion is powered only by the core stage. Differently from the VLS-1, this launch vehicle does not perform a non-powered coast-arc. Figure 8 presents the altitude and relative velocity profiles.

Figure 8
(a) Altitude profile of the Ariane 5 and (b) velocity profile of the Ariane.

From Fig. 8 we can notice that the behavior of the altitude profile firstly exceeded the desired altitude (LEO with 200 km) reaching a maximum altitude, and then the desired altitude is obtained. The reason for that is presumably because the flight does not perform a non-powered phase, so the vehicle takes longer to get the right inclination in order to injected into orbit.

VLS-ALFA LAUNCH VEHICLE

It is known that the VLS-Alfa will replace the last two solid stages of the former VLS-1 by a single liquid upper stage. Then since the VLS-Alfa is an improvement of the former VLS-1, the VLS-1 will be used as a reference vehicle to perform the simulations. The upper stage of the VLS-Alfa presumably will perform a coast phase, so the L75 is supposed to support restart capability. The mission is to launch a satellite into a reference circular orbit of 500 km of altitude from the Alcântara Launch Center (2°22'39.52" S, 44°23'57.71" W). The parameters of the vehicles are given in Tables 5 and 6.

Table 5
Key parameters: Brazilian launch vehicle VLS (Schlingloff 1987).

VLS	mprop (kg)	ms (kg)	Isp (s)	tb (s)	CD*S (m2)
1st Stage	28900	6200	257.90	62.826	3.0
2nd Stage	7140	1680	279.10	62.087	1.3
3rd Stage	4370	1330	270.74	58.267	0.0
4th Stage	820	170	281.85	74.546	0.0

Table 6
Key parameters: Brazilian launch vehicle VLS-Alfa.

VLS-Alfa	mprop (kg)	ms (kg)	Isp (s)	tb (s)	CD*S (m2)
1st Stage	28900	6200	257.90	62.826	3.0
2nd Stage	7140	1680	279.10	62.087	1.3
3rd Stage (before coasting)	4370	1022.53	315.00	243.657	0.0
4th Stage (after coasting)	820	1022.53	315.00	46.219	0.0

As the mission of the VLS-Alfa is still not totally defined, the propellant mass of the upper stage had to be estimated. Thus, an amount of 6900 kg was estimated based on internal reports. From this value an amount of 1100 kg was taken for the phase after coasting, i.e., for orbit injection. The altitude and relative velocity profiles for both vehicles are presented in Fig. 9.

Figure 9
(a) Relative velocity profile of the former and future Brazilian launch vehicle and (b) altitude profile of the former and future Brazilian launch vehicle.

CONCLUSIONS

This article presented a modeling of a launch vehicle using object-oriented programming. Two models for trajectory optimization were presented. To allow a better visualization of the codes, the Unified Modeling Language (UML) was used. The future Brazilian VLS-Alfa and its predecessor VLS-1 were simulated and compared. Both of them are supposed to support perform coast phase. The performance of the Ariane 5, which does not perform coast phase, was presented and the results showed good agreement with the literature. Since the stages of the launch vehicle is an assembly engine and propellant system, the modular approach is interesting to allow future integration of the propulsion system in order to study the coupling among the disciplines. Accordingly, the modular approach using object-oriented programming can facilitate a multidisciplinary optimization.

REFERENCES

• Balesdent M (2011) Multidisciplinary design optimization of launch vehicles (PhD Dissertation). Nantes: Ecole Centrale de Nantes.
• Betts JT (1998) Survey of numerical methods for trajectory optimization. Journal of Guidance, Control, and Dynamics 21(2):193-207. doi: 10.2514/2.4231
  » https://doi.org/10.2514/2.4231
• Brauer GL, Cornick DE, Stevenson R (1977) Capabilities and applications of the program to optimize simulated trajectories (POST). (CR-2770). NASA Contractor Report.
• Brown KR, Harrold EF, Johnson GW (1969) Rapid optimization of multiple-burn rocket flights. (CR-1430). NASA Contractor Report.
• Cornelisse JW, Schöyer HFR, Wakker KF (1979) Rocket propulsion and spaceflight dynamics. London: Pitman.
• Gath PF, Calise AJ (2001) Optimization of launch vehicle ascent trajectories with path constraints and coast arcs. Journal of Guidance, Control and Dynamics 24(2):296-304. doi: 10.2514/2.4712
  » https://doi.org/10.2514/2.4712
• Hargraves CR, Paris SW (1987) Direct trajectory optimization using nonlinear programming and collocation. Journal of Guidance, Control and Dynamics 10(4):338-342. doi: 10.2514/3.20223
  » https://doi.org/10.2514/3.20223
• Herman AL, Conway BA (1996) Direct optimization using collocation based on high-order Gauss-Lobatto quadrature rules. Journal of Guidance, Control and Dynamics 19(3):592-599. doi: 10.2514/3.21662
  » https://doi.org/10.2514/3.21662
• Hinckel JN (1995) An object oriented approach to launch vehicle performance analysis. Presented at: 31st AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit; San Diego, USA. doi: 10.2514/6.1995-3094
  » https://doi.org/10.2514/6.1995-3094
• Jacob HG (1972) An engineering optimization method with application to STOL-aircraft approach and landing trajectories. (TN-D-6978). NASA Technical Report.
• Miele A, Wang T (2003) Multiple-subarc gradient-restoration algorithm, part 2: application to a multistage launch vehicle design. Journal of Optimization Theory and Applications 116(1):19-39. doi: 10.1023/A:1022154001343
  » https://doi.org/10.1023/A:1022154001343
• Pontani M, Teofilatto P (2014) Simple method for performance evaluation of multistage rockets. ACTA Astronautica 94(1):434-445. doi: 10.1016/j.actaastro.2013.01.013
  » https://doi.org/10.1016/j.actaastro.2013.01.013
• Rao AV (2009) A survey of numerical methods for optimal control. Advances in Astronautical Sciences 135(1):497-528.
• Schlingloff H (2005) Astronautical engineering: an introduction to the technology of spaceflight. Bad Abach: Ingenieurbuero Dr. Schlingloff Publications.
• Schlingloff H (1991) SKYNAV a design tool for space launcher analysis and flight performance optimization. Munich: User's Documentation.
• Schlingloff H (1987) Control laws for optimal spacecraft navigation. Journal of Spacecraft and Rockets 24(1):48-51. doi: 10.2514/3.25871
  » https://doi.org/10.2514/3.25871
• Seywald H (1994) Trajectory Optimization Based on Differential Inclusion. Journal of Guidance, Control, and Dynamics 17(3):480-487. doi: 10.2514/3.21224
  » https://doi.org/10.2514/3.21224
• Silva CSC (1995) Simulação de Sistemas de Motores Foguetes a Propelentes Líquidos (Master's Dissertation). São José dos Campos: Instituto Tecnológico de Aeronáutica. In Portuguese.
• Teren F, Spurlock OF (1966) Payload optimization of multistage launch vehicles. (TN D-3191). NASA Technical Note.
• Tewari A (2007) Atmospheric and spaceflight dynamics. Boston: Birkhäuser.
• von Stryk O, Bulirsch R (1992) Direct and indirect methods for trajectory optimization. Annals of Operations Research 37(1):357-373. doi: 10.1007/BF02071065
  » https://doi.org/10.1007/BF02071065

Edited by

Section Editor: Luiz Martins-Filho

Publication Dates

• Publication in this collection
  13 Sept 2018
• Date of issue
  2018

History

• Received
  08 June 2017
• Accepted
  06 Nov 2017