*Article* **Analysis on the Isostatic Bipod Mounts for the HERA Mission LIDAR**

### **Nicole G. Dias 1, Paulo Gordo 1,2, Hugo Onderwater 2, Rui Melicio 3,4,\* and António Amorim <sup>1</sup>**


**Abstract:** The Light Detection and Ranging (LIDAR) is a time-of-flight altimeter instrument being developed for the HERA mission, designated as Planetary ALTimeter (PALT). PALT is positioned in the center of the top face of the HERA probe, and therefore, it cannot use radiators to stabilize its internal temperature. The contribution of this paper is the design of isostatic bipod mounts for the LIDAR primary mirror. The performance of PALT must be maintained over a wide operational range, from −60 ◦C to 80 ◦C. These temperature requirements imply that a careful isostatic mount structure design is critical to maintaining performance in all operational scenarios. The purpose of the instrument is to perform range measurements from 500 m to 14 km. The instrument will contribute to the detailed characterization of the asteroid's topography, assist the probe navigation in operations such as fly-bys (including on the dark side of the asteroid) or landing. PALT has an emitter system that generates 2 ns, 100 μJ, 1535 nm laser pulses and a receiver system that collects the backscattered signal from the asteroid. The receiver system is composed of a 70 mm diameter Cassegrain telescope and a refractive system that focuses the signal on the sensor.

**Keywords:** isostatic bipod mounts; LIDAR; PALT; HERA mission; a-thermalization; optical performance

### **1. Introduction**

The HERA spacecraft includes several payload instruments, such as the Time-of-Flight (ToF) LIDAR that will measure the distances from the HERA spacecraft to the target. The measurement operations shall be performed at a distance from 500 m to 14 km, enabling operations such as fly-bys or landings. Previous space missions have deployed analogous instruments for specific requirements. One of the main challenges in those missions was the operational temperature range, since the LIDAR instruments were directly exposed to space and the required optical tolerances to maintain the instrument performance, i.e., internal alignment of optics and alignment between receiver and emitter. Each mission requires a specific LIDAR measurement range, operational temperature interval, radiation requirements and target objects, making the LIDAR design rather specific. The contribution of this paper is the design of isostatic bipod mounts for a small mirror of the LIDAR.

Within the mirror, the reflecting coating, the substrate and the supports have to be validated over the thermal range of the mission. When the mirror is subjected to a thermal condition, maintaining the surface shape is critical. "A temperature change induces a change in size which, for curving optical surfaces, produces a change in radius and hence a focus shift" [1]; therefore, a large change in temperature can cause a distortion.

Isostatic mounts, or flexures, are optomechanical components capable of compensating thermal and gravitational deformations in optical components. Isostatic bipod mounts are

**Citation:** Dias, N.G.; Gordo, P.; Onderwater, H.; Melicio, R.; Amorim, A. Analysis on the Isostatic Bipod Mounts for the HERA Mission LIDAR. *Appl. Sci.* **2022**, *12*, 3497. https://doi.org/10.3390/ app12073497

Academic Editors: Simone Battistini, Filippo Graziani and Mauro Pontani

Received: 28 January 2022 Accepted: 27 March 2022 Published: 30 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

used to prevent deformations due to differences of thermal expansion coefficients, and to keep the mirror's optical axis in place [1]. Several design limitations occur due to the mission environmental requirements.

An a-thermal design must be taken into consideration since critical components have a different Coefficient of Thermal Expansion (CTE). Usually, optical components have a low CTE compared to the mount, thus the a-thermal design is a critical aspect to the mount design [2]. Pijnenburg [2] studied two different a-thermalization methods: one matching the bond thickness, and the second by applying elastic elements in the mount. The first approach is to adapt the thickness of the adhesive to compensate for the expansion of the optics and mounts. However, it resulted in the instability of the optical components, since it resulted in a thick bond line which does not work properly for strength and stiffness. The second approach included a tangential isostatic flexible element, a leaf spring. If well dimensioned, the mirror would be supported by three leaf springs, with a small bonding spot. The stability under inertial loads, such as gravity, is achieved by constraining with high stiffness the degree of freedom of the optical component, where "only residual local stresses in the glass arises around the bond spots" [2].

Following the second approach of Pijnenburg [2], several designs were presented for several space missions, especially for mirrors with large aperture. There are two types of conformities that should be followed: radial compliance for axisymmetric mirrors to compensate for thermal expansion mismatch, and tangential compliance to prevent assembly stress from navigating towards the mirror surface. Lateral mirror supports on the edge of the mirror help in fulfilling these compliances [3,4]. Isostatic mounts can be categorized according to the type of flexure element: simple blade flexures, that are used for tangential edge support for small axisymmetric mirrors, and a bipod flexure, which is a combination of two blade flexures forming a triangle [3,4]. Kihm [3] presented a design for a lightweight primary mirror, with pockets in the back surface, and with three square bosses extruded at the edge of the mirror for isostatic mounting. The isostatic mount presented by Kihm [3] had three different components: the flexure A, which is permanently bonded onto the mirror's boss; the flexure B, which is fastened to flexure A and is the support flexure; and shims that are placed between the flexures. In [5], an adjustable bipod flexure for a large aperture mirror is presented, formed by two flexure bars with tangential and radial blades. The connection between the flexure and the mirror is made by an invar connector since the invar's CTE matches with the material of the mirror, and it is glued to the mirror with epoxy adhesive.

In this paper, the isostatic bipod mounts design of PALT is presented. It follows a similar approach to the flexure presented by Kihm [3,4], which means that the apex of the triangle formed by the flexure should point to the center of mass of the mirror to minimize surface distortion [3,4]. Each bipod flexure has a symmetric combination of two tangential blades and the radial blade: tangential blades give tangential compliance, and radial blades give radial compliance to the mirror. Thus, when the mirror is accelerated in lateral directions or perpendicular to the mirror's optical axis, the tangential flexures compensate for this effect; however, radial expansion of the mirror due to thermal loads can be compensated for by the radial flexures [4]. The major differences from previous implementations lie in the mirror size and in the isostatic mount bipod mechanical design.

This paper is structured as follows: Section 2 presents the previous research. Section 3 presents the PALT design. Section 4 explains the assembly and bonding procedures. Section 5 reveals the thermoelastic simulations and respective results. Finally, in Section 6, the conclusions are outlined.

#### **2. Previous Research**

The Planetary ALTimeter (PALT) mechanical design has been optimized on different aspects to meet the requirements established for the HERA mission. In [6–8], the HELENA LIDAR's performance was evaluated regarding weaker requirement of shock, static simulation of 5G in two different directions, and thermal conductive simulations between −40 ◦C

and +60 ◦C. The main objective was to study the secondary mirror deformation towards the base of the telescope, and maximum stresses. Two different materials (aluminum and titanium) were assessed for the upper part of the LIDAR telescope. The simulations revealed a higher factor of safety, and a minor displacement of the secondary mirror towards the base of the telescope for the titanium material. This modification was implemented on PALT.

In order to guarantee that the best optical performance is maintained, the primary mirror shall be a-thermalized, which is possible by the optimizing the isostatic bipod mounts (see Figure 1). Furthermore, the assembly and alignment procedures of the primary mirror should be as simple as possible to avoid damages on the mirror. This paper presents an optimization of the isostatic bipod mounts design, and details on the assembly procedure.

**Figure 1.** (**a**) HELENA design; (**b**) HELENA engineering model. (**c**) PALT design. Several transformations were made since HELENA design, with the isostatic bipod mounts being a critical design improvement.

### **3. Design**

The design was iterated several times, with the focus being on the initial design and a final design. The first design of the isostatic bipod mounts was presented in HELENA [6,7], which was the starting point of the isostatic bipod mount for PALT. This first design is depicted in Figure 2.

**Figure 2.** First design of isostatic bipod mount, which are positioned symmetrically around the mirror, and the representation of the adhesive pads area (highlighted in green).

This design followed the one implemented on HELENA LIDAR [6,7] (see Figure 3).

**Figure 3.** Isostatic bipod mount assembled on HELENA.

The primary mirror of HELENA had a 100 mm diameter and 11 mm thickness.

For PALT, a 70 mm diameter mirror was required. A first approach on the isostatic bipods mount design was based on the HELENA's optomechanics. The simulations of this design are presented in Section 4, and it was concluded that further optimization was desirable.

In a second version, the mirror's thickness was decreased to 8 mm and the three extrude bosses were equally distributed with an 8 mm length (where isostatic bipod mounts are assembled). Both isostatic bipod mounts combine tangential and radial blades in each support leg, arranged in a triangle. The proposed design includes two blade radial flexures and a tangential blade between the radial blades, whose purpose is to compensate for radial and tangential deformations, respectively. The material chosen for this isostatic mount is Titanium Grade 5 (Ti-6Al-4V) (see Table 1).


**Table 1.** Properties of the materials.

The bipod isostatic mount presented in this article is shown in Figure 4, and it has dimensions of 33 mm height and a maximum length of 8 mm. Since the optimized implementation is a key element to achieve the required performance, the optimization of the design of the isostatic bipod mount was divided into four elements: a base flexure, two columns, and a top fitting. The separation of the bipod eases the bonding process between the mounts and the mirror extrude bosses.

**Figure 4.** Isostatic bipod mount that will be assembled in PALT. They are positioned symmetrically around the mirror, and the adhesive pads area are displayed (highlighted in blue).

The extruded areas presented in the inner upper structure of the isostatic bipod mount are adhesive areas where the mirror is to be attached. The adhesive is presented as a thin pad with 0.1 mm thickness.

### **4. Assembly and Bonding Procedures**

In HELENA, the isostatic bipod mount was bonded to the mirror with a 100 μm layer of adhesive. The bonding process was implemented by using a 100 μm diameter nylon wire, between the isostatic mount and the mirror, which were removed after the adhesive was fully cured. It was verified that the mount design blocked the application of the adhesive, which sprung the need for design changes. The isostatic bipod mount updated design is easily assembled and, since it is divided into four individual pieces, allow a better application of the epoxy without compromising the mirror. The bonding procedure is still to be qualified; however, the bonding procedure and the alignment procedure is defined. The first step fixates the mirror in space using three supports, assisted by Kapton pads and adjustment screws. A sheet made of Kapton, is to be placed on the mirror extrusions to ensure the adhesive is kept on the designated area and not touching the mirror optics surface (see Figure 5).

**Figure 5.** Mirror fixation and protection. The protective sheet is perforated to allow an easy removal after assembly and bonding of the flexures.

The adhesive thickness remains 100 μm apart from the flexures, on all four sides. To accomplish this, Kapton tape of a specific thickness is to be placed on the corners of the mirror extrusions. The tape serves, not only to guarantee that all four sides of the mirror extrusion will remain 100 μm apart from the flexures, but also to prevent any excess adhesive that may spread from the beyond the designated area (highlighted surfaces depicted in Figure 4). All the protective bonding components will be removed before the end of the work life of the adhesive.

After assuring that the mirror is not contaminated by adhesive, the flexure base is placed under the mirror extrusion, with the adhesive already applied. The columns of the flexure are placed on the flexure base, both with the adhesive already applied, finalizing with the top lid of the flexure, with the adhesive already applied, and screwed into place, fixating the upper part of the flexure (see Figure 6). The two screws guarantee the connection of the isostatic bipod mount: the screw is inserted on the top fitting, goes through the columns, and reaches the flexure base.

**Figure 6.** Final step on the bonding procedure. The protective sheets are removed before the end of the work life of the adhesive (schematic).

### **5. Thermoelastic Simulations**

This section presents the thermoelastic simulations of both designs and the respective results.

### *5.1. Thermal Modeling and Materials*

The isostatic bipod mounts were simulated in simple conditions to assess their behavior. The simulation only contemplated a circular base, to maintain the system symmetric, the three isostatic bipod mounts, and the mirror. In Figure 7, the simulation performed for the first design is depicted, and in Figure 8, the simulation performed for the optimized flexure design. The boundary conditions applied to the model were to prevent axial and tangential displacement, which means the base only has freedom radially, and a fixed-point support on the base center, which is coincident with the optical axis. The temperature was defined in the bottom surface of the base.

**Figure 7.** The boundary conditions of the first design.

**Figure 8.** The boundary conditions of the final design.

Both designs were subjected to the qualification non-operational temperature limits, i.e., +80 ◦C and −60 ◦C. The materials chosen were also the same: titanium grade 5 (Ti-6Al-4V) for the isostatic bipod mounts, aluminum 7075-T7351 for the base, Zerodur for the mirror, and Epoxy 2216B/A Gray for the adhesive pads. The properties of the chosen materials are presented in Table 1. The mechanical screw behavior was not assessed at this stage.

The Poisson's ratio and the Young's Modulus were given by the supplier and are confidential. The remaining properties are based on documentation [9,10]. The properties presented in Table 1 are the inputs needed for the thermo-elastic simulation.

The tolerance of the system is given by an optical study (i.e., tolerance analysis in ZEMAX). Regarding the telescope, it was identified that the worst tolerance offenders are the mirror radius and the mirror displacement. The tolerances that are assessed in this paper is the displacement on XZ plane of the primary mirror, the distance of the primary mirror to the base, and the curvature radius of the primary mirror (see Table 2).


**Table 2.** Thermal acceptable tolerances.

#### *5.2. Thermoelastic Results*

The aim of these simulations was to assess the isostatic bipod mount's behavior in the extreme temperatures. Two topics were studied, namely the survivability of the mirror (i.e., assessment of the stress in the mirror) and the analysis of the mirror displacement regarding its nominal position (to evaluate if the flexures were optimized and working correctly).

Figures 9 and 10 represent the XZ displacement results for the first isostatic bipod mount design for the hot and cold non-operational scenarios, respectively.

**Figure 9.** Result of hot (+80 ◦C) non-operational scenario—First design: (**a**) *X*-axis displacement, (**b**) *Z*-axis displacement.

**Figure 10.** Result of cold (−60 ◦C) non-operational scenario—First design: (**a**) *X*-axis displacement, (**b**) *Z*-axis displacement.

From the results, it is evident that there is a displacement of the mirror relative to its nominal position. From the hot scenario, the extreme values of XZ displacement are in range of −3.461 × <sup>10</sup><sup>−</sup>4; 3.027 × <sup>10</sup>−<sup>4</sup> mm, and for the cold scenario, the extreme values of XZ displacement are in range of −4.280 × <sup>10</sup><sup>−</sup>4, 4.894 × <sup>10</sup>−<sup>4</sup> mm, which is fully within the range of acceptable tolerances from Table 2. An animation of the simulations can be found in (https://www.youtube.com/watch?v=M6bY21nINoE accessed on 23 March 2020)

Regarding the stress, the critical components were analyzed, namely the mirror and the isostatic bipod mounts. The thermal conductance along the mirror and the adhesives, and between the adhesives and the isostatic bipod mounts is established as 2500 W/(m2K). The contacts between the other materials (aluminum and titanium) are established as 150 W/(m2K) (as defined in [11]). The results of the von-Mises stress for the hot non-operational case and for the cold non-operational case are presented in Figures 11 and 12, respectively.

**Figure 11.** Result of hot (+80 ◦C) non-operational scenario—First design: (**a**) Mirror stress, (**b**) Isostatic bipod mount stress (deformation scale factor: 76).

**Figure 12.** Result of cold (−60 ◦C) non-operational scenario—First design: (**a**) Mirror stress, (**b**) Isostatic bipod mount stress (deformation scale factor: 76).

The mirror maximum stress is sited on the bonding areas, with a safety factor of 3.29 for the hot scenario, and 2.33 for the cold scenario, which means that the mirror will withstand the non-operational temperatures.

The isostatic bipod mounts have a good overall performance, since their displacement is located on the radial blades. The tangential blade does not present a considerable stress since the displacement condition is essentially radial. The maximum stress is located on the base of the isostatic bipod mounts, but since the critical stress is on the blades and in the bonding area, that maximum is not considered. When the maximum stress range is decreased, the value of stress starts to appear on the radial blade, which is coincident with the radial deformation of the base (see Figures 11b and 12b). The safety factor for the isostatic bipod mount considers the maximum stress value, which converts into a safety factor of 2.82 for the hot scenario, and 1.99 for the cold scenario.

Making a similar study for the final design of the isostatic bipod mounts, the XZ displacement results are depicted in Figures 13 and 14 for the hot and cold non-operational scenarios, respectively.

**Figure 13.** Result of hot (+80 ◦C) non-operational scenario—Final design: (**a**) *X*-axis displacement, (**b**) *Z*-axis displacement.

From the results, it is evident that there is a displacement of the mirror relative to its nominal position. Nevertheless, the results are within the established values for XZ displacement tolerances from Table 2. The extreme values of the hot non-operational temperatures for the final design of the isostatic bipod mounts are within the range of −1.159 × <sup>10</sup><sup>−</sup>4; 1.133 × <sup>10</sup>−<sup>4</sup> mm, and for the cold non-operational temperatures are in range of −1.602 × <sup>10</sup><sup>−</sup>4; 1.683 × <sup>10</sup>−<sup>3</sup> mm.

Regarding the stress in the mirror and the isostatic bipod mounts, the same assessment previously mentioned was made. The results are presented in Figures 15 and 16 for the hot and cold non-operational scenario, respectively.

**Figure 15.** Result of hot (+80 ◦C) non-operational scenario—Final design: (**a**) Mirror stress, (**b**) Isostatic bipod mount stress (deformation scale factor: 73).

The mirror's maximum stress is located on the bonding areas, with a safety factor of 1.75 for the hot scenario, and 1.24 for the cold scenario. Since the safety factor of the mirror is higher than unity, the mirror can withstand the hot and cold non-operational scenarios.

Regarding the isostatic bipod mounts stress distribution, it can be concluded that they have a good overall performance, since their displacement is located on the radial blades. The tangential blade does not present a considerable stress since the displacement condition is essentially radial.

The safety factor for the isostatic bipod mount converts into a safety factor of 3.66 for the hot scenario, and 2.59 for the cold scenario. Since the safety factor is higher than unity, the isostatic bipod mounts can withstand the hot and cold non-operational scenarios.

The maximum stress is found on the designated area of the radial blades (see Figures 15b and 16b), which is expected.

The distance between the primary mirror and the base is given in Figures 17 and 18.

**Figure 18.** Result of *Y*-axis displacement for the final design: (**a**) hot (+80 ◦C) non-operational scenario, (**b**) cold (−60 ◦C) non-operational scenario.

The distance of the mirror towards the base for the first design is between [−0.0248, 0.035] mm, and for the final design is between [−0.0340, 0.0480] mm. These values are within the requirement of ±0.09 mm defined in Table 2.

For the analysis of curvature radius variation, the following result processing approach was implemented: (1) the displacement offset of the mirror was removed; (2) the difference between the mirror outer edge and mirror inner edge was calculated (see Figure 19). To this difference we named the radius curvature variation of the mirror, and its results are presented in Table 3.

**Figure 19.** Representation of the inner and outer edge of the mirror: (**a**) first design, (**b**) final design. **Table 3.** Curvature radius deformation.


According to Table 2, a 1 mm curvature radius corresponds to a surface displacement of ±2 μm. The cold scenario does not fulfil this requirement; and for this case, the LIDAR performance has a loss of energy.

### **6. Conclusions**

The paper describes the development of an isostatic bipod mount design for small mirrors that withstand severe temperature requirements. The assessment of the design was presented considering a preliminary design, which was already built but not yet tested, and the updated design, which is going to be built and subjected to testing.

The first design presents several complicated processes for bonding and alignment, which prompt a design that would be easier to integrate. Both designs present a viable design for the LIDAR, having taken into account the XZ displacement, the distance of the mirror towards the base, the curvature radius, and the stress on the mirror. Both designs survive the non-operational scenarios; however, the final design provides an easier implementation and bonding procedure, since it is divided into four individual parts, in contrast to the first design. Additionally, it is expected that the final design has a higher resistance to shock, because the bonding of the mirror is implemented in four perpendicular areas per isostatic bipod mount instead of only one.

The future steps include the assessment of the performance of the isostatic bipod mounts assembled in the PALT, and consequently, the evaluation of the optical performance of PALT.

**Author Contributions:** Conceptualization, N.G.D., P.G., H.O., R.M. and A.A.; methodology, N.G.D., P.G., H.O., R.M. and A.A.; software, N.G.D.; investigation, N.G.D., P.G., H.O., R.M. and A.A.; writing original draft preparation, N.G.D.; writing—review and editing, N.G.D., P.G., H.O., R.M. and A.A.; validation, N.G.D., P.G., H.O., R.M. and A.A.; supervision, N.G.D., P.G., H.O., R.M. and A.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** European Union's Horizon 2020 research and innovation programme under Grant Agreement No. 870377 (Project No. NEO-MAPP). FCT (Foundation for Science and Technology, I.P.), through: CENTRA, project UIDB/00099/2020; LAETA project UIDB/50022/2020; ICT project UIDP/04683/2020.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Nomenclature**


### **References**


**Paolo Teofilatto 1,\*, Stefano Carletta <sup>1</sup> and Mauro Pontani <sup>2</sup>**


**\*** Correspondence: paolo.teofilatto@uniroma1.it

**Abstract:** This papers introduces an analytic method to define multistage launcher trajectories to determine the payload mass that can be inserted in orbits of different semimajor axes and inclinations. This method can evaluate the gravity loss, which is the main term to be subtracted to the Tziolkowski evaluation of the velocity provided by the thrust of a launcher. In the method, the trajectories are dependent on two parameters only: the final flight-path angle *γ<sup>f</sup>* at the end of the gravity-turn arc of the launcher trajectory and the duration *tc* of the coasting arc following the gravity-turn phase. The analytic formulas for the gravity-turn phase, being solutions of differential equations with a singularity, allow us to identify the trajectory with a required final flight-path angle *γ<sup>f</sup>* in infinite solutions with the same initial vertical launch condition. This can also drive the selection of the parameters of the pitch manoeuvre needed to turn the launcher from the initial vertical arc. For any pair *γ<sup>f</sup>* and *tc*, a launcher trajectory is determined. A numerical solver is used to identify the values *γ<sup>f</sup>* and *tc*, allowing for the insertion of the payload mass into the required orbit. The analytic method is compared with a numerical code including the drag effect, which is the only effect overlooked in the analytic formulas. The analytical method is proven to predict the payload mass with an error never exceeding the 10% of the actual payload mass, found through numerical propagation.

**Keywords:** multistage launch vehicles; ascent trajectory optimization; analytical performance evaluation; rocket staging

### **1. Introduction**

The number of microsatellites to be orbited is increasing rapidly, in particular for the delivery and replacement of microsatellites as part of mega constellation programs. This is stimulating a new wave of dedicated launch vehicles capable of offering responsive, flexible and cost-effective services to this huge and new market. Some of these "micro-launchers" are listed in Table 1. The mass of payload *mpay* to be inserted into LEO is below 1000 kg. In fact, *mpay* is the driving element of a space rocket and it is important to define it at the first step of a launcher design. A very popular method of computation of *mpay* is offered by the Tziolkowski formula. However, the results obtained using this method are not accurate since the space environment (gravity and aerdoynamic forces) and the guidance (direction of thrust) are not taken into account. On the other hand, the accurate computation of a launcher trajectory and the evaluation of launcher performances is a typical and complex problem in the optimization of aerospace trajectories. Launcher performances are evaluated maximizing the final mass that can be set into orbits with different parameters (semimajor axis, eccentricity, inclination). The problem has initial- and final-state variable constraints, path constraints and discontinuities in the mass variation, due to stage separation. This optimization problem can be approached in different ways: all the approaches require educated guesses to initialize some iterative algorithm able to refine the initial guesses and to achieve the optimal solution.

**Citation:** Teofilatto, P.; Carletta, S.; Pontani, M. Analytic Derivation of Ascent Trajectories and Performance of Launch Vehicles. *Appl. Sci.* **2022**, *12*, 5685. https://doi.org/10.3390/ app12115685

Academic Editor: Jérôme Morio

Received: 21 April 2022 Accepted: 30 May 2022 Published: 3 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).


**Table 1.** Listof Micro-Launchers.

In the scientific literature, some papers [1–7] address the problem of optimizing the ascent trajectory of launch vehicles through indirect approaches. Examples are the multiple-subarc gradient restoration algorithm, proposed by Miele, and multiple shooting techniques [5–7]. The previously cited works [2–7] require a considerable deal of effort for deriving the analytical conditions for optimality and for the subsequent programming and debugging. Furthermore, these methods can suffer from a slow rate of convergence and an excessive dependence on the starting guess. Alternative approaches [8,9] are direct in nature, often are more robust, but require discretization of the problem (e.g., through collocation [8,9]) and the use of dedicated algorithms for nonlinear programming, without avoiding the need of a starting guess. Most recently, the indirect heuristic method [10] was proposed as an effective approach for trajectory optimization of ascent vehicles, with the remarkable feature of not requiring any starting guess. For the purpose of preliminary analysis of the performance of new ascent vehicles, fast, approximate approaches that are exempt from convergence issues and minimize the computational effort are desirable. With this intent, Pontani and Teofilatto [11] and Pallone et al. [12,13] propose two numerical approaches for performance evaluation of multistage launch vehicles that aim at obtaining accurate predictions of the performance attainable from a launch vehicle.

The aim of this paper is to derive formulas, possibly as simple as the Tziolkowski one, able to obtain in a few seconds the performance of a multistage space rockets and also able to offer trajectory parameters useful as input for the convergence of complex numerical programs with accurate representation of the rocket system and its space environment.

The formulas originate from previous work [14], which has been generalized to multistage rockets in [15]; all the rocket state variables are derived in [16]. In the present paper, the method is extended to include the optimal choice of a cost arc and the guidance of the last stage of the launcher till the injection into the selected orbit, so the performance of a space launcher is obtained. The method is based on simple mathematics and physics principles and only a few seconds are needed to determine the payload mass that can be set into orbits of different altitude and inclinations, starting from any location of the launch pad.

The major simplification made in order to obtain these results is that the aerodynamic drag is neglected, and this produces an error in the computation of the launcher trajectories. It is rather difficult to treat the aerodynamic effect analytically due to its highly non-linear dependence on velocity and altitude. Some analytic formulas including aerodynamics are derived in [17]; these however are restricted to the vertical arc of the trajectory. To evaluate the impact of drag, the method presented was compared to numerical results obtained by a high-fidelity level of representation of the launcher dynamics. It is shown that the error in the estimate of the optimal payload mass to be inserted into orbit is relevant only for orbits

of altitude below 600 km, and in these cases the error does not exceed 10% Moreover, it is proved that the analytic approach proposed here can provide good guesses for numerical solvers of the optimal control problem.

The paper is organized as follows: in Section 2, a first sizing of a space launcher is given based on the Tziolkowski formula. Assuming a payload mass to be delivered in LEO orbit and assuming a fixed percentage of losses, a method is derived to size the staging of a multistage rocket in an optimal way. In Section 3, analytical formulas for the launcher gravity-turn phase are developed in detail for the case of multiple stage rockets. These formulas provide the launcher trajectory of the first stages and allow for the evaluation of the so-called gravity loss. Sections 4 and 5 describe the method of computation of the full launcher trajectory till the injection into orbit. Section 6 is dedicated to testing the method with respect to numerical computations.

### **2. First Sizing of a Launch Vehicle**

Achieving a single stage to orbit is still an infeasible task and space launchers have a number *N* of stages ranging from 2 to 4. Each stage has a total mass *mi*

$$m\_i = m\_{p\_i} + m\_{si}$$

where *mpi* is the propellant contained within the tank of each stage and the structural mass coefficient of each stage is defined as

$$
\varepsilon\_{\bar{l}} = \frac{m\_{s\_{\bar{l}}}}{m\_{\bar{i}}} \tag{1}
$$

The total mass at lift off is

$$M\_1 = \Sigma\_{i=1}{}^N m\_i + \left(m\_{\text{lss}} + m\_{\text{pay}}\right)^\*$$

with *mhs* the mass of the heat shields (here released at the end of the gravity-turn phase) and *mpay* is the satellite mass to be injected into orbit. At the time *tb*<sup>1</sup> of the end of the first stage boost (burn time), the propellant mass *mp*<sup>1</sup> is consumed and the mass of the space rocket is

$$M\_{1f} = M\_1 - m\_{p\_1}$$

The velocity reached at burn time *tb*<sup>1</sup> can be computed integrating the acceleration

$$\dot{V} = \frac{T}{m}\cos\alpha - \frac{\mu}{r^2}\sin\gamma - \frac{1}{2}\rho V^2 \frac{\mathbb{C}\_D}{m} \tag{2}$$

here *T* is the thrust intensity with a direction determined by the angle *α* taken from the velocity *V*- . The velocity direction is determined by the flight-path angle *γ* taken from the local horizon direction ˆ *θ*, see Figure 1. The variable *r* is the radius distance from the centre of the Earth to the launcher centre of mass and *μ* is the Earth gravitational constant. The atmospheric density is denoted by *ρ* and *CD* is the launcher drag coefficient computed with respect to the reference surface *S*.

To underline the thrust's contribution to the variation in velocity and consider all the other terms as "losses", let us write (2) as

$$\dot{V} = \frac{T}{m} - \frac{\mu}{r^2} \sin \gamma - \frac{1}{2} \rho V^2 \frac{C\_D S}{m} - \frac{T}{m} \left(1 - \cos \alpha \right)^2$$

Integrating the above acceleration between time zero and *tb*1, one obtains the velocity after the first-stage boost

$$
\Delta V\_1 = \int\_{t\_0}^{t\_{b1}} \frac{T}{m} \, dt - \int\_{t\_0}^{t\_{b1}} \frac{\mu}{r^2} \sin \gamma \, dt - \int\_{t\_0}^{t\_{b1}} \frac{1}{2} \rho V^2 \frac{\mathbb{C}\_D S}{m} \, dt - \int\_{t\_0}^{t\_{b1}} \frac{T}{m} \left(1 - \cos \alpha \right) \, dt \tag{3}
$$

**Figure 1.** The state variable of the launch trajectory. The inertial frame (*r*ˆ0, ˆ *θ*0) has origin on the centre of the Earth (in dark blue). The orbit frame (*r*ˆ, ˆ *θ*) has origin on the launcher centre of mass. The red arrow defines the velocity and the blue arrow the thrust direction

Equation (3) is called the "equation of losses"

$$
\Delta V\_1 = \Delta V\_{Tz} - \Delta V\_{\text{graw}} - \Delta V\_{drag} - \Delta V\_{\text{wis}}.
$$

where the first variation of velocity due to the engine boost is called the Tziolkowski variation of velocity and the other velocities subtracted to Δ*VTz* are called gravity, drag and misalignment losses, respectively. Average thrust is defined by *T* = *gIsp m*˙ , where *m*˙ is equal to the mass rate *<sup>m</sup>*˙ <sup>=</sup> <sup>−</sup>*dm dt* . Inserting *<sup>T</sup>* <sup>=</sup> <sup>−</sup>*gIsp dm dt* in (3), we obtain the Tziolkowski velocity variation provided by the first-stage engine

$$
\Delta V\_{Tz} = gIsp \, \log(\frac{M\_1}{M\_{1f}})
$$

Introducing the subrocket structural mass ratios

$$M\_k = \frac{M\_{kf}}{M\_k},\ \ M\_k = \Sigma\_{i=k}{}^N m\_i + \left(m\_{\text{par}} + m\_{\text{hs}}\right),\ \ M\_{kf} = M\_k - m\_{p\_k},\ \ k = 1, N$$

the variation of the velocity provided by the engines of a multistage rocket is

$$
\Delta V\_{Tz} = \Sigma\_{k=1}^{N} \otimes Isp\_k \log\left(\frac{1}{\mathcal{U}\_k}\right) \tag{4}
$$

The values of the specific impulses *Ispk* depend on the kind of propellant, for instance liquid or solid. Figure 2 shows the specific impulses for stages with solid propellant: these values are taken from the database [18] and we have average values

$$Isp\_1 = 280\,\text{s} \,, \, Isp\_2 = 290\,\text{s} \,, \, Isp\_3 = 270\,\text{s} \tag{5}$$

Other parameters defined by technical constraints are the structural mass ratio of each stage (1): Figure 3 shows these values for a number of stages taken from the database [18], and we have average values

$$
\epsilon\_1 = 0.11, \; \epsilon\_2 = 0.13, \; \epsilon\_3 = 0.3 \tag{6}
$$

**Figure 2.** The specific impulse of different stages taken from the database ©Mark Wade [18].

**Figure 3.** The structural mass ratio of different stages taken from the database ©Mark Wade [18].

It is important to know that if the total mass *M*1, the payload mass *mpay* and the values *k*ith the five unknow, *Ispk* are given, there is an optimal selection of *Uk*, maximizing the Tziolkowski velocity (4). The procedure to identify the optimal subrocket mass ratios *U*∗ *<sup>k</sup>* is in Appendix A.

The parameters *k*, *Ispk* can be easily guessed at the first step of the rocket design, and the maximum payload mass *mpay* to be inserted into a specific reference orbit is the performance characterizing a space launcher.

With *k*, *Ispk*, and *mpay* fixed, the total mass *M*<sup>1</sup> determines the optimal mass distribution *U*∗ *<sup>k</sup>* (staging) and the maximum Tziolkowski velocity, according to Appendix A. If this velocity is equal to the velocity required to inject the payload mass *mpay* into the reference orbit, *M*<sup>1</sup> is the launcher lift-off mass achieving this performance.

One point to underline is that the velocity to be provided by the thrust is not equal to the orbital velocity of a spacecraft on a circular orbit of altitude *h*: *Va* = *μ <sup>h</sup>*+*RE* . In fact, we must consider the losses: these depend on each phase of the ascent trajectory. In the present Section it is assumed that the launcher velocity due to thrust must be overcome by 30% the required value for orbit injection. The orbit velocity relative to the Earth *V*-*R* depends on the latitude *L* of the launch site and on the azimuth *ψ*

$$V\_R = \sqrt{V\_a^2 + V\_E^2 - 2V\_a V\_E \sin \psi}$$

where *VE* = *ω<sup>E</sup> RE* cos(*L*) is the velocity of the launch site and the azimuth angle *ψ* is related to the launch site latitude and orbit inclination by the equation: sin *<sup>ψ</sup>* <sup>=</sup> cos *<sup>i</sup>* cos *<sup>L</sup>* .

Then, the required velocity is here defined as

$$V\_{req} = V\_R(1 + 30\%)$$

Figure 4 shows the relative (in blue) and required (in red) velocities to achieve a reference circular orbit of altitude h = 650 km with different inclinations from a launch site at L = 30 deg N. Table 2 shows the minimum mass at lift off for a micro-launcher having parameters (5) and (6) for different payload masses (from 100 to 1000 kg) to be injected into the reference orbit (*Vreq* = 9805 m/s). The ratio *<sup>M</sup>*<sup>1</sup> *mpay* is constant.

**Table 2.** Minimumlift-off mass to inject the payload mass into circular polar orbit of altitude *h* = 650 km (average parameters (5) and (6), and launch site L = 30 deg).


**Figure 4.** Therelative (dotted) and required (solid line) velocities to achieve circular orbits of altitude h = 800 km (black), 700 km (green), 600 km (blue), 500 km (red) and different inclinations

### **3. A Second Step (Gravity Losses)**

At lift off, the launcher performs a vertical trajectory followed by a pitch manoeuvre to select the required orbit plane. Then, to limit the aerodynamic loads on the structure, any launcher is forced to follow a zero angle of attack trajectory (*α* = 0) during the atmospheric flight. This trajectory is called the gravity-turn trajectory, and it is performed by the first stages of space launchers. The trajectory is basically planar and the state variables can be the relative velocity, flight-path angle, altitude and range (*V*, *γ*, *h*,*s*). The equations of motion are

$$\begin{cases} \dot{V} = \frac{T}{m} - \frac{\mu}{r^2} \sin \gamma - \frac{1}{2} \rho V^2 \frac{C\_D S}{m} \\ \dot{\gamma} = \frac{\dot{s}}{R\_E} - \frac{\mu}{r^2} \cos \gamma + \frac{1}{2} \rho V \frac{C\_L S}{m} \\ \dot{h} = V \sin \gamma \\ \dot{s} = \frac{R\_E}{r} V \cos \gamma \end{cases} \tag{7}$$

In fact, these equations are singular under the lift-off condition; *<sup>V</sup>* = 0, *<sup>γ</sup>*<sup>0</sup> = *pi* 2 .

As observed in [14], this singularity generates an infinite number of solutions corresponding to the same lift-off condition: one can use this singularity to select among the different solutions the one with a prescribed final value, for instance a fixed flight-path angle *γ<sup>f</sup>* [15].

An analytic formula for the velocity *V*, due to Culler and Fried, was obtained under the following hypotheses:


The last of these is the strongest assumption, and is handled considering an average value *n* of the thrust-to-weight ratio. This average value can be computed using the formulas for the variation of mass *m* = *m*<sup>0</sup> − *md t* = *m*0(1 − *q t*) and the definition of the mass ratio *U* = *ms m*<sup>0</sup> , which verifies *U* = 1 − *q tb*.

Then

$$\begin{aligned} \bar{m} &= \quad \frac{1}{t\_b} \int\_0^{tb} \frac{T}{m\mathcal{G}}(t) \, dt = \frac{T}{\mathcal{G}} \int\_0^{tb} \frac{dt}{1 - qt} = \\ &= \quad \frac{n\_0}{t\_b} \frac{1}{-q} \left[ \ln(1 - qt) \right]\_0^{t\_b} = \frac{n\_0}{-q} \ln(1 - qt\_b) = 0 \\ &= \quad -\frac{n\_0}{1 - U} \ln(U) = \frac{n\_0}{1 - U} \ln\left(\frac{1}{U}\right) \end{aligned}$$

where *n*<sup>0</sup> = *<sup>T</sup> gm*<sup>0</sup> is the initial thrust-to-weight ratio.

Now two transformations are defined:


By the transformation in Equation (1) we have the following differential equations for the velocity and kick angle

$$\begin{cases} \dot{V} = -\underline{g}\cos\chi + \underline{g}\overline{n} \\ \dot{\chi} = \stackrel{\circ}{\theta}\sin\chi \end{cases} \tag{8}$$

Now the following trigonometrical equations are recalled

$$\begin{aligned} \cos \chi &= \begin{array}{c} 1 - \tan^2 \frac{\chi}{2} \\ 1 + \tan^2 \frac{\chi}{2} \end{array} = \frac{1 - z^2}{1 + z^2} \\\ \sin \chi &= \sin \left( 2 \frac{\chi}{2} \right) = 2 \sin \frac{\chi}{2} \cos \frac{\chi}{2} = 2 \tan \frac{\chi}{2} \cos^2 \frac{\chi}{2} = \frac{2z}{1 + z^2} \end{array}$$

showing that the transformation to the *z*-variable allows us to transform the trigonometric function in rational functions. Then, Equation (8) becomes

$$\begin{cases} \dot{V} = -g\,\vec{n} - g\frac{1-z^2}{1+z^2} \\ \dot{z} = \frac{g}{V}z \end{cases} \tag{9}$$

Dividing the first by the second equation of the system (9), one has

$$
\frac{dV}{dz} = \bar{n}\frac{V}{z} - \frac{V}{z}\frac{1-z^2}{1+z^2}
$$

$$
\frac{dV}{V} = \left(\frac{\bar{n}}{z} - \frac{1}{z}\frac{1-z^2}{1+z^2}\right)dz
$$

that is

By integration

$$V = A \, z^{n-1} (1 + z^2) \tag{10}$$

The velocity in (10) is a function of the variable *z*. The time dependence can be obtained deriving the relation between time and *z*. Substitution of (10) in the second equation of the system (9) gives

$$\dot{z} = \frac{dz}{dt} = -\frac{\mathcal{S}}{A} \frac{1}{z^{n-2}(1+z^2)}$$

by separation of variables and integrating from *t*<sup>0</sup> = 0 to *t*

$$t = \frac{A}{\mathcal{g}} \left( \frac{z^{n-1}}{\bar{n} - 1} + \frac{z^{n+1}}{\bar{n} + 1} \right) \tag{11}$$

The launch initial conditions are in the given hypothesis

$$\begin{cases} V\_0 = 0\\ \gamma\_0 = \frac{\pi}{2} \\\ \chi\_0 = 0 \Rightarrow z\_0 = 0 \end{cases}$$

Imposing these initial conditions in the solution (10), one has

$$V\_0 = A \, z\_0^{n-1} (1 + z\_0^2) = 0 \quad \forall A$$

Then, the initial conditions are verified for any value of the constant *A*. That is, there are infinite (gravity-turn) solutions for the velocity, parameterized by A, corresponding to the same initial condition. This violation of the Chauchy theorem on uniqueness of the solution of a differential equation is due to the singularity of *z*˙ for *V* = 0 in the system (9).

One can take advantage of this multiplicity of solutions to select for instance the solution that achieves a chosen flight-path angle *γ<sup>f</sup>* at the burn-out time *tb*. The chosen *γ<sup>f</sup>* defines a value *<sup>z</sup> <sup>f</sup>* = tan( *<sup>π</sup>*/2−*γ<sup>f</sup>* <sup>2</sup> ). By Equation (11), one obtains

$$A = \lg t\_b \left( \frac{z\_f^{n-1}}{\bar{n} - 1} + \frac{z\_f^{n+1}}{\bar{n} + 1} \right)^{-1} \tag{12}$$

and the related velocity is obtained inserting this value of *A* into Equation (10).

In particular, the velocity at burn out is

$$V\_f = A \, z\_f^{n-1} (1 + z\_f^2) \tag{13}$$

In addition, the altitude *h* and range *s* can be easily derived [16]. For the altitude

$$\dot{h} = V \sin \gamma = V \cos \chi = V \frac{1 - z^2}{1 + z^2}$$

Then

$$\frac{dh}{dz} = \frac{h}{\dot{z}} = \frac{V^2}{gz} \frac{1 - z^2}{1 + z^2}$$

By Equation (10), one has

$$\frac{d\mathbf{h}}{dz} = \frac{A^2}{\mathcal{S}} \left( z^{2n-3} - z^{2n+1} \right).$$

The elementary integration of the above equation between current *z* and initial *z*<sup>0</sup> gives the formula for the altitude

$$h(z) = h(z\_0) + \frac{A^2}{g} \left[ \frac{z^{2n-2}}{2\overline{n} - 2} - \frac{z^{2n+2}}{\overline{n} + 2} \right]\_{z\_0}^{z\_f} \tag{14}$$

For the range *s*

$$\text{ ș } = V \cos \gamma = V \sin \chi = 2\, V \frac{z}{1+z^2}$$

then

$$\frac{ds}{dz} = \frac{\dot{s}}{\dot{z}} = 2\frac{A^2}{\mathcal{S}} \left( z^{2\mathfrak{n}-1} + z^{2\mathfrak{n}+1} \right)$$

The integration gives the formula for the range

$$s(z) = s(z\_0) + 2\frac{A^2}{g} \left[ \frac{z^{2n}}{2\bar{n}} + \frac{z^{2n+2}}{2\bar{n}+2} \right]\_{z\_0}^{z\_f} \tag{15}$$

Equations (10), (12), (14) and (15), together with (11), provide the state variables of the gravity-turn trajectory of a single-stage rocket with flight-path angle *γ<sup>f</sup>* at burn-out time *tb*.

For instance, let us consider initial lift off condition (*V*<sup>0</sup> = 0, *χ*<sup>0</sup> = 0, *h*<sup>0</sup> = 0,*s*<sup>0</sup> = 0) and choose the final condition *γ<sup>f</sup>* = 0 this condition corresponds to *z <sup>f</sup>* = 1. Then Equation (12) gives

$$A = \lg t\_b \frac{\vec{n}^2 - 1}{2\vec{n}}$$

By Equations (10), (14) and (15), the velocity, altitude and range at burn out are equal to

$$V\_f = 2A = \lg t\_b \frac{\bar{n}^2 - 1}{\bar{n}} \tag{16}$$

$$h\_f = \frac{A^2}{\mathcal{g}} \frac{1}{\bar{n} - 1} = \frac{\mathcal{g}t\_b^2}{2} \frac{\bar{n}^2 - 1}{\bar{n}^2}$$

$$s\_f = \frac{A^2}{\mathcal{g}} \left(\frac{1}{\bar{n}} + \frac{1}{\bar{n} + 1}\right) = \frac{\mathcal{g}t\_b^2}{2} \left(\frac{\bar{n}^2 - 1}{\bar{n}}\right)^2 \left[\frac{1}{2\bar{n}} + \frac{1}{2\bar{n} + 2}\right]$$

Note that the final velocity formula of a rocket is generally obtained by the Tziolkowski formula that takes into account only the thrust acceleration. In the present setting, the Tziolkowski formula gives the final velocity

$$V\_{Tz} = \mathcal{g} \,\,\tilde{m}t\_b\,\,$$

The final velocity obtained here, which includes gravity effect, is given by (16): the difference among the two is the gravity loss and it is equal to

$$
\Delta V\_{\mathcal{S}} = V\_{Tz} - V\_f = \frac{\mathcal{g} \, t\_b}{\vec{n}} \tag{17}
$$

The above results can be generalized to multistage launchers.

### **4. Analytic Derivation of Launcher Trajectories**

The above analytic formulas can be generalized to multistage launchers. Let us start with two stages with the following average parameters (*n*1, *Isp*1, *u*1) and (*n*2, *Isp*2, *u*2).

The first stage has variables *z*1, *V*<sup>1</sup> = *A*<sup>1</sup> *z n*¯ <sup>1</sup>−1 <sup>1</sup> (<sup>1</sup> + *<sup>z</sup>*<sup>2</sup> <sup>1</sup>) in the time interval [*t*<sup>10</sup> = 0, *t*<sup>1</sup> *<sup>f</sup>* = *tb*<sup>1</sup> ]. The second stage has variables *z*2, *V*<sup>2</sup> = *A*<sup>2</sup> *z n*¯ <sup>2</sup>−1 <sup>2</sup> (<sup>1</sup> + *<sup>z</sup>*<sup>2</sup> <sup>2</sup>) in the time interval [*t*<sup>20</sup> = *tb*<sup>1</sup> , *t*<sup>2</sup> *<sup>f</sup>* = *tb*<sup>2</sup> + *tb*<sup>1</sup> ] and the final value *z*2*<sup>f</sup>* is fixed (see Figure 5).

Of course the two solutions*z*1, *V*<sup>1</sup> and *z*2, *V*<sup>2</sup> must match at the boundary point

$$z\_{1f} = z\_{20}$$

$$V\_{1f} = V\_1(z\_{1f}) = V\_2(z\_{20}) = V\_{20}$$

The two matching conditions imply

$$\frac{A\_1}{A\_2} = z\_{20}^{n\_2 - n\_1} \tag{18}$$

Equation (11) gives at first-stage burn out

$$t\_{b1} = \frac{A\_1}{g} \left( \frac{z\_{20}^{n\_1 - 1}}{\bar{n}\_1 - 1} + \frac{z\_{20}^{n\_1 + 1}}{\bar{n}\_1 + 1} \right) = ^{def} f\_1(z\_{20}) \tag{19}$$

At the burn out of the second stage, one has *z*<sup>2</sup> *<sup>f</sup>* and

$$t\_{b2} = \frac{A\_2}{g} \left( \frac{z\_{2f}^{\eta\_2 - 1}}{\bar{n}\_2 - 1} + \frac{z\_{2f}^{\eta\_2 + 1}}{\bar{n}\_1 + 1} - \frac{z\_{20}^{\bar{n}\_2 - 1}}{\bar{n}\_2 - 1} + \frac{z\_{20}^{\bar{n}\_2 + 1}}{\bar{n}\_1 + 1} \right) = ^{def} f\_2(z\_{20}) \tag{20}$$

We have a system of three Equations (18)–(20) in the three unknowns *A*1, *A*<sup>2</sup> and *z*20. To solve it, let us consider the ratio

$$\frac{t\_{b2}}{t\_{b1}} = \frac{f\_1}{f\_2}$$

Because of (18), this ratio is a function of *z*<sup>20</sup> only

$$f(z\_{20}) = \frac{t\_{b2}}{t\_{b1}} - \frac{f\_1}{f\_2} = 0$$

The solution *z*<sup>20</sup> of *f*(*z*20) = 0 (in the range [0,1]) provides the value of the flightpath angle at the end of the first stage. Moreover, by Equation (20), one finds *A*<sup>2</sup> and by Equation (19) one finds *A*1. Hence, the velocity profile of the two-stage launcher with a prescribed flight-path angle *z*<sup>2</sup> *<sup>f</sup>* at second-stage burn out is equal to

$$\begin{cases} V(z) = A\_1 z\_1^{\eta\_1 - 1} (1 + z\_1^2) \ for \ z \in [0, z\_{20}] \\ V(z) = A\_2 z\_2^{\eta\_2 - 1} (1 + z\_2^2) \ for \ z \in [z\_{20}, z\_{2f}] \end{cases} \tag{21}$$

With same arguments, we obtain the altitude and range

$$\begin{aligned} h(z) &= \frac{A\_1^2}{\mathcal{S}} \left[ \frac{z\_1^{2n\_1 - 2}}{2\bar{n}\_1 - 2} - \frac{z\_1^{2n\_1 + 2}}{\bar{n}\_1 + 2} \right], \text{ for } z \in [0, z\_{20}]\\ h(z) &= h(z\_{20}) + \frac{A\_2^2}{\mathcal{S}} \left[ \frac{z\_2^{2n\_2 - 2}}{2\bar{n}\_2 - 2} - \frac{z\_2^{2n\_2 + 2}}{\bar{n}\_2 + 2} \right]\_{z\_{20}}^z, \text{ for } z \in [z\_{20}, z\_{2f}] \end{aligned} \tag{22}$$

$$\begin{aligned} s(z) &= 2\frac{A\_1^2}{g} \left[ \frac{z^{2n\_1}}{2\bar{n}\_1} + \frac{z^{2n\_1+2}}{2\bar{n}\_1+2} \right], \text{ for } z \in [0, z\_{20}]\\ s(z) &= R(z\_{20}) + 2\frac{A\_2^2}{g} \left[ \frac{z^{2n\_2}}{2\bar{n}\_2} + \frac{z^{2n\_2+2}}{2\bar{n}\_2+2} \right]\_{z\_{20}}^z, \text{ for } z \in [z\_{20}, z\_{2f}] \end{aligned} \tag{23}$$

Note that the gravity loss of the two-stage launcher is

$$
\Delta V\_{\mathcal{S}} = gIsp\_1 \log(\frac{1}{\mathcal{U}\_1}) + gIsp\_2 \log(\frac{1}{\mathcal{U}\_2}) - V(z\_{2f}),
$$

Any number of stage N can be handled: of course, more numerical work is needed to find out the unknown variables. The N-1-matching conditions on the flight-path angle generate the following N-1-matching conditions on the velocity

$$\begin{array}{rcl} V\_1(z\_{20\prime}A\_1) &=& V\_2(z\_{20\prime}A\_2) \\ V\_2(z\_{30\prime}A\_2) &=& V\_3(z\_{30\prime}A\_3) \\ &\vdots \\ V\_{N-2}(z\_{N-2\cdot 0\prime}A\_{N-2}) &=& V\_{N-1}(z\_{N-1\cdot 0\prime}A\_{N-1}) \end{array}$$

These, together with the *N* relationships related to each burn time

$$\begin{array}{rcl} t\_{b\_1} &=& t\_{b\_1}(A\_1, z\_{20}) \\ t\_{b\_2} &=& t\_{b\_2}(A\_2, z\_{20}, z\_{30}) \\ &\vdots \\ t\_{b\_N} &=& t\_{b\_N}(A\_N, z\_{20}, \dots, z\_{N0}) \end{array}$$

generate 2*N* − 1 equations to be solved with respect to the 2*N* − 1 unknowns *z*20, *z*30, ... *zN* 0, *A*1,... *AN*.

By these solutions, the velocity, altitude and range profiles can be computed using Equations (13)–(15) in the ranges [0, *z*20],... [*zN* 0, *z <sup>f</sup>* ].

For instance, for a three-stage launcher, one has

$$t\_{b1} = \frac{A\_1}{\mathcal{S}} \left( \frac{z\_{20}^{n\_1 - 1}}{n\_1 - 1} + \frac{z\_{20}^{n\_1 + 1}}{n\_1 + 1} \right) \tag{24}$$

$$t\_{b2} = \frac{A\_2}{g} \left[ \left( \frac{z\_{30}^{n\_2 - 1}}{n\_2 - 1} + \frac{z\_{30}^{n\_2 + 1}}{n\_2 + 1} \right) - \left( \frac{z\_{20}^{n\_2 - 1}}{n\_2 - 1} + \frac{z\_{20}^{n\_2 + 1}}{n\_2 + 1} \right) \right] \tag{25}$$

$$t\_{b3} = \frac{A\_3}{\mathcal{S}} \left[ \left( \frac{z\_3^{n\_3 - 1}}{n\_3 - 1} + \frac{z\_3^{n\_3 + 1}}{n\_3 + 1} \right) \right]\_{z\_{\mathcal{B}\_0}}^{z\_{\mathcal{B}\_f}} \tag{26}$$

$$A\_1 = A\_2 z\_{20}^{\overline{n}\_2 - \overline{n}\_1} \tag{27}$$

$$A\_2 = A\_3 \ z\_{30}^{\overline{n}\_3 - \overline{n}\_2} \tag{28}$$

(these five Equations (24)–(28) with the five unknowns *A*1, *A*2, *A*3, *z*20, *z*30).

One solution is to write the variables *A*<sup>1</sup> and *A*<sup>2</sup> as function of *A*<sup>3</sup> using Equations (27) and (28). Then, the ratios *tb*<sup>2</sup> *tb*1 , *tb*3 *tb*1 are obtained dividing (25) and (26) by the Equation (24). The two ratios give two equations on the two unknown *z*20, *z*30. The solutions can be found numerically within the interval [0,1], and they determine the values of the flight-path angles at the end of the first and second stages. The values *A*1, *A*2, *A*<sup>3</sup> are derived by Equations (24) and (26) and and the velocities at the end of each stage are obtained.

For a four-stage launcher, there are seven equations in the seven unknowns *A*1, *A*2, *A*3, *A*4, *z*20, *z*30, *z*<sup>40</sup>

$$A\_1 = A\_2 z\_{20}^{\overline{\pi}\_2 - \overline{\pi}\_1} \tag{29}$$

$$A\_2 = A\_3 \, z\_{30}^{\eta\_3 - \eta\_2} \tag{30}$$

$$A\_3 = A\_4 \ z\_{40}^{\mathfrak{m}\_4 - \mathfrak{m}\_3} \tag{31}$$

$$t\_{b1} = \frac{A\_1}{\mathcal{g}} \left( \frac{z\_{20}^{n\_1 - 1}}{n\_1 - 1} + \frac{z\_{20}^{n\_1 + 1}}{n\_1 + 1} \right) \tag{32}$$

$$t\_{b2} = \frac{A\_2}{\mathcal{S}} \left[ \left( \frac{z\_{30}^{n\_2 - 1}}{n\_2 - 1} + \frac{z\_{30}^{n\_2 + 1}}{n\_2 + 1} \right) - \left( \frac{z\_{20}^{n\_2 - 1}}{n\_2 - 1} + \frac{z\_{20}^{n\_2 + 1}}{n\_2 + 1} \right) \right] \tag{33}$$

$$t\_{b\_3} = \frac{A\_3}{\mathcal{S}} \left[ \left( \frac{z\_{40}^{n\_3 - 1}}{n\_3 - 1} + \frac{z\_{40}^{n\_3 + 1}}{n\_3 + 1} \right) - \left( \frac{z\_{30}^{n\_3 - 1}}{n\_3 - 1} + \frac{z\_{30}^{n\_3 + 1}}{n\_3 + 1} \right) \right] \tag{34}$$

$$t\_{b\_4} = \frac{A\_4}{g} \left[ \frac{2}{n\_4^2 - 1} - \left( \frac{z\_{40}^{n\_4 - 1}}{n\_4 - 1} + \frac{z\_{40}^{n\_4 + 1}}{n\_4 + 1} \right) \right] \tag{35}$$

The first three equations allow us to write *A*1, *A*2, *A*<sup>3</sup> as function of *A*<sup>4</sup> only. Then, the three ratios *tb*<sup>2</sup> *tb*1 , *tb*3 *tb*1 , *tb*4 *tb*1 are functions of the three unknowns *z*20, *z*30, *z*40. The solution is used then to derive *A*<sup>1</sup> − *A*4.

Note that the above formulas for three- or four-stage launchers provide the values of the flight-path angles and the velocities at the end of each stage for the launcher trajectory with a fixed value of the final flight-path angle.

#### **5. Guidance**

The gravity-turn trajectory constrains the guidance and should be abandoned if possible; the relevant condition is the factor *q α*, where *q* = <sup>1</sup> <sup>2</sup> *<sup>ρ</sup> <sup>V</sup>*<sup>2</sup> is the dynamical pressure and *α* is the angle of attack. This factor must be below a limit [*q α*]*max* which is peculiar to any launcher and can be satisfied by the gravity-turn (*α* = 0) trajectory. As the altitude increases, the atmospheric density *ρ* decreases exponentially, so at a certain point of the flight the condition of *q α* disappears and a different guidance can be applied. The guidance applied here is based on: (A) a possible coasting arc, and (B) thrust along the horizontal direction.

Generally, a coasting arc is performed after the burn out of the stage *N* −1 of an *N* stage launcher to increase the altitude at the cost of a moderate decrease in velocity. After the coasting arc, the launcher is close to the target altitude *hf* , whereas the required velocity *Vreq* must be gained. This velocity basically has the horizontal direction ˆ *θ* = cos *γ V*ˆ + sin *γ* ˆ *l*, then the thrust direction is along ˆ *θ* during the N-stage boost.

With the above guidance scheme, the full launcher trajectory depends on two parameters: the flight-path angle *γ<sup>f</sup>* at the end of the gravity-turn phase and the duration *tc* of the coasting arc. The two parameters (*γ<sup>f</sup>* , *tc*) are introduced into an iterative routine (such as the Matlab routine fsolve.m) and updated to set the errors on the final altitude and velocity (*h*(*tb*) − *hf* , *V*(*tb*) − *Vreq*) to zero, where *tb* is the N-stage burn-out time.

Let us consider a three-stage rocket as an example. The first two stages are in gravity turn with a prescribed final flight-path angle *γ*<sup>2</sup> *<sup>f</sup>* , then a costing arc of duration *tc* is performed. Finally the third-stage boost is applied with horizontal thrust direction. The Matlab routine fsolve.m is used to update the values (*γ*<sup>2</sup> *<sup>f</sup>* , *tc*), for instance, to achieve the polar circular orbit of altitude *hf* = 650 km from a launch site with a latitude of 30 deg N. The required velocity is then equal to *VR* = 7542 m/s.

The launcher has the parameters reported in Table 3. The Tziolkowski velocity is equal to Δ*VTz* = 9743 m/s. From the data of Table 3 it is possible to derive the mass parameters of each stage *ms*<sup>1</sup> , *mp*<sup>1</sup> , *ms*<sup>2</sup> , *mp*<sup>2</sup> , *ms*<sup>3</sup> , *mp*<sup>3</sup> and the mass of each substage *M*2, *M*3, see Equations (A1) and (A2) in Appendix A. From the values *n*01 , *n*<sup>02</sup> , *n*<sup>03</sup> and *M*1, *M*2, *M*3, it is possible to derive the propellant mass rates *m*˙ 1, *m*˙ 2, *m*˙ 3, hence the burn times *tbi* <sup>=</sup> *mpi <sup>m</sup>*˙ *<sup>i</sup>* , *<sup>i</sup>* = 1, 3. For any fixed *<sup>γ</sup>*2*<sup>f</sup>* , the Equations (21)–(23) give the state of the launcher at the end of the second stage where the gravity-turn phase ends. After a coasting of duration *tc*, the third stages thrust the rocket up to the burn time *tb*<sup>3</sup> ; the final altitude and velocity are compared with the required (*hf* , *VR*). The parameters *γ*2*<sup>f</sup>* , *tc* are updated in order to match the required final orbit conditions. After few iterations, the routine finds the following parameters to reach the required orbit: *γ*<sup>2</sup> *<sup>f</sup>* = 26.5 deg, *tc* = 340 s.

**Table 3.** Dataof the three-stage launcher.


### **6. Numerical Comparison**

The results obtained by the analytic formulas are now compared with a numerical code simulating the launcher trajectory with better accuracy. The equations of motion are the three-dimensional equations of flight in the relative velocity frame with state variables (*r*, *λ*, *L*, *VR*, *γR*, *ψR*)

$$\begin{aligned} \dot{r} &= V\_R \sin \gamma\_R\\ \dot{\lambda} &= V\_R \cos \gamma\_R \frac{\cos \psi\_R}{r \cos \theta} L\\ \dot{L} &= \frac{V\_R \cos \gamma\_R \sin \psi\_R}{r} \\ \dot{V}\_R &= r \omega\_E^2 \cos L(\cos L \sin \gamma\_R - \sin L \cos \gamma\_R \sin \psi\_R) + \frac{f\_V}{m} \\ &\dot{\gamma}\_R = \frac{V\_R}{r} \cos \gamma\_R + 2\omega\_E \cos L \cos \psi\_R + \\ &+ r \frac{\omega\_E^2}{V\_R} \cos L(\cos L \cos \gamma\_R + \sin L \sin \gamma\_R \sin \psi\_R) + \frac{f\_I}{mV} \\ \dot{\psi}\_R &= -\frac{V\_R}{r} \tan L \cos \gamma\_R \cos \psi\_R + 2\omega\_E \cos L \tan \gamma\_R \sin \psi\_R + \\ &- r \frac{\omega\_E^2}{V\_R \cos \gamma\_R} \cos L \sin L \cos \psi\_R - 2\omega\_E \sin L \end{aligned} \tag{36}$$

with

$$\begin{aligned} \frac{f\_V}{m} &= \frac{T}{m} \cos \alpha - \frac{\mu}{r^2} \sin \gamma\_R - \frac{1}{2} \rho \, \mathcal{C}\_D \frac{S}{m} V\_R^2\\ \frac{f\_I}{m} &= \frac{T}{mV} \sin \alpha - \frac{\mu}{V\_R r^2} \cos \gamma\_R + \frac{1}{2} \rho \, \mathcal{C}\_L \frac{S}{m} V\_R \end{aligned} \tag{37}$$

Equations (36) and (37) are used with *α* = 0 for the gravity-turn arc, with *T* = 0 for the coasting arc and with *α* = *γ* for the third-stage boost. Lift off begins with an initial vertical arc of small duration *tV*; this is computed in the Cartesian frame with non-singular coordinates. After the vertical arc, to turn the trajectory from the vertical direction, a pitch manoeuvre starts with a fixed direction of thrust *α* = *θ<sup>k</sup>* and ends when *γ<sup>R</sup>* = *θk*. At this time (called pitch over), the gravity-turn trajectory starts. The pitch direction of thrust *θ<sup>k</sup>* is chosen so to have the prescribed value *γ*2*<sup>f</sup>* at the second-stage burn out. Then, the coasting arc duration is chosen to achieve the required orbit.

Figures 6–8 show the analytical and numerical graphs of the velocity, flight-path angle and altitude during gravity turn. Note that at the end of the second stage the discrepancy between analytic and numeric results is rather small for the velocity, whereas a reduction of 15% occurs in the altitude due to the numerical integration of the drag force. In the numerical calculation the atmospheric density is represented by the USATMO 76 model and the drag coefficient is kept constant at an average value *CD* = 0.4. The reference surface is *S* = *<sup>π</sup> <sup>d</sup>*<sup>2</sup> <sup>4</sup> , with a diameter of *d* = 2 m. The drag effect reduces the payload mass that can be inserted into the reference orbit.

**Figure 6.** The velocity plot: analytic and numeric results during the gravity-turn arc. Numerical velocity of the first and second stage is in black, analytic velocity of the first stage is in blue and second stage in orange line.

**Figure 7.** The flight-path angle plot: analytic and numeric results during gravity-turn arc. Numerical flight-path angle of the first and second stage is in black, analytic flight-path angle of the first stage is in blue and second stage in orange line.

In fact, the numerical results give insertion in the polar circular orbit of altitude 650 km of a payload of mass *mpay* = 390 kg, vertical flight duration *tv* = 5 s, and pitch angle *θ<sup>k</sup>* = 76 deg (generating *γ*2*<sup>f</sup>* = 26.5 deg), and coasting time *tc* = 355 s. The payload mass reduction with respect to the analytical results is about 4%. The evolution of the orbit parameters during the launch is shown in Figures 9 and 10, and Figure 11 shows the launcher ground track.

**Figure 8.** The altitude plot: analytical and numerical results during gravity-turn arc. Numerical altitude of the first and second stage is in black, analytical altitude of the first stage is in the blue line and the second stage in the orange line.

**Figure 9.** Thein-plane orbit parameters during the ascent. First-stage parameters are the blue line; coasting and second stage are the orange line.

**Figure 10.** Out–of-plane orbit parameters during the ascent. First-stage parameters are in the blue line; coasting and second stage are in the orange line.

**Figure 11.** Launcher ground track.

Figure 12 shows the payload mass evaluation for polar circular orbits of different altitudes obtained by the analytic (blue curve) and numerical (orange curve) algorithms. Of course, the bigger differences are for orbits of lower altitude; however, the reduction in the payload mass does not exceed 10% of the value found analytically.

**Figure 12.** Payload mass of polar circular orbits of different altitudes, analytical (blue line) vs numeric results (orange line).

### **7. Concluding Remarks**

This research extends the preceding generalization of the Culler and Fried formulation to ascent trajectories that include a coast arc and the guidance of the upper stage. The analytical developments related to the multistage formulation of the Culler and Fried analysis are contained in Section 4, dedicated to closed-form expressions for some fundamental parameters of multistage launch vehicles. Section 5 is focused on inclusion of the coast arc and guidance of the upper stage. These steps lead to the defining of an analytical methodology for the preliminary evaluation of the performance and trajectory of a launch vehicle, whose preliminary mass distribution is optimized through the solution of a polynomial equation of the same order as the number of stages (cf. Appendix A). The comparison of the analytical method at hand with the ascent path obtained through numerical integration testifies its accuracy, particularly in terms of prediction of the final payload mass, especially

for medium- and high-altitude final circular orbits, with the error never exceeding the 10% of the actual payload mass.

**Author Contributions:** Conceptualization, P.T. and M.P.; methodology, P.T.; software, P.T.; validation, S.C. and M.P.; formal analysis, P.T.; investigation, P.T., S.C. and M.P.; resources, P.T., S.C. and M.P.; data curation, P.T., S.C. and M.P.; writing—original draft preparation, P.T.; writing—review and editing, S.C. and M.P.; supervision, P.T.; project administration, P.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Appendix A**

Let us recall the definitions of stage and subrocket structural mass ratios of an N-stage launcher:

$$\epsilon\_k = \frac{m\_{s\_k}}{m\_k} \quad \mathcal{U}\_k = \frac{M\_{k - m\_{p\_k}}}{M\_k} \quad k = 1, \dots, N$$

where *mk* = *msk* + *mpk* is the mass of each stage (structure plus propellant mass) and *Mk* = Σ*<sup>N</sup> <sup>j</sup>*=*kmk* + *mpay* is the mass of each subrocket.

$$m\_{s\_k} = \frac{M\_k(1 - \mathcal{U}\_k)}{1 - \epsilon\_k}, \; m\_{p\_k} = M\_k(1 - \mathcal{U}\_k) \tag{A1}$$

$$M\_k = M\_{k-1} \frac{\left(\mathcal{U}\_{k-1} - \epsilon\_{k-1}\right)}{1 - \epsilon\_{k-1}} \ , k = 2 , N - 1 \ , \ M\_N = m\_{\text{pay}} \tag{A2}$$

Let the values *M*<sup>1</sup> and *mpay* be given, then, by

$$\frac{m\_{\text{pay}}}{M\_1} = \Pi\_{i=2}^N \frac{M\_i}{M\_{i-1}} \tag{A3}$$

The Tziolkowski velocity is equal to

$$
\Delta V\_{Tz} = -\lg \sum\_{i=1}^{N} Isp\_i \, \log(lI\_i),
$$

The problem is now: select the values *U*∗ *<sup>k</sup>* to minimize

$$J = \operatorname{g} \Sigma\_{i=1}^{N} \operatorname{Isp}\_{i} \log(\mathcal{U}\_{i})$$

under the constraint (A3), which is equal to (see (A2)):

$$\Phi = \frac{m\_{\text{ray}}}{M\_1} - \Pi\_{i=2}^N \frac{\mathcal{U}\_i - \mathfrak{e}\_i}{1 - \mathfrak{e}\_i} = 0 \tag{A4}$$

To solve the constrained minimum problem, the Hamilton function is introduced with Lagrangian multiplier *λ*:

$$H = f + \lambda \,\Phi$$

and the solution for the *N* + 1 unknowns *Ui*, *λ* derives from the *N* + 1 equations

.

$$\begin{aligned} \frac{\partial H}{\partial M\_i} &= 0\\ \frac{\partial H}{\partial \lambda} &= 0 \end{aligned} \tag{A5}$$

These equations correspond to

$$\frac{gIsp\_i}{lI\_i} - \frac{\lambda}{1 - \varepsilon\_i} \Pi\_{k=1, k \neq i}^{N} \frac{lI\_k - \varepsilon\_k}{1 - \varepsilon\_k} = 0\tag{A6}$$

and to Equation (A4). From (A6), the Lagrangian multiplier *λ* is derived:

$$\lambda = \frac{g \operatorname{Isp}\_i(\mathcal{U}\_i - \mathfrak{e}\_i)}{\mathcal{U}\_i} \Pi\_{k=1}^N \frac{1 - \mathfrak{e}\_k}{\mathcal{U}\_k - \mathfrak{e}\_k} \tag{A7}$$

Equation (A7) holds true for any *i* = 1, *N*, and this corresponds to

$$\frac{g\,I\,sp\_1(\mathcal{U}\_1-\epsilon\_1)}{\mathcal{U}\_1} = \dots \dots = \frac{g\,I\,sp\_N(\mathcal{U}\_N-\epsilon\_N)}{\mathcal{U}\_N}$$

that is

$$\frac{\text{g}\,\text{Isp}\_{i}(\text{}\,\text{l}\,\text{l}\_{i}-\text{e}\_{i})}{\text{l}\,\text{l}\_{i}} = \frac{\text{g}\,\text{Isp}\_{1}(\text{}\,\text{l}\,\text{l}\_{1}-\text{e}\_{1})}{\text{l}\,\text{l}\_{1}}\,\text{,}\,\text{i}=\text{2,N}\tag{\text{A8}}$$

Equation (A8) corresponds to

$$\mathcal{U}I\_{i} = \frac{Isp\_{i}\,\varepsilon\_{i}}{\mathcal{U}\_{1}(Isp\_{i} - Isp\_{1}) + Isp\_{1}\,\varepsilon\_{1}} \mathcal{U}\_{1} \,\,\, i = 2, N \tag{A9}$$

Replacing the Formulas (A9) into (A4), one has the following algebraic equation of order N in the unique unknown *U*1:

$$\frac{m\_{\rm pay}}{M\_1} \Pi\_{i=1}^N \left[ \mathcal{U}\_1 (Isp\_i - Isp\_1) + Isp\_1 \epsilon\_1 \right] - \Pi\_{i\_1}^N \frac{\epsilon\_i}{1 - \epsilon\_i} (Isp\_1 (\mathcal{U}\_1 - \epsilon\_1))^N = 0 \tag{A10}$$

The solution *U*∗ <sup>1</sup> of (A10) is used in (A9) to obtain all the subrocket structural mass ratios \$ *U*∗ *i* % *<sup>i</sup>*=1,*N*, maximizing the Tziolkowski variation of velocity Δ*VTz*.

### **References**

