Emission-Driven Hybrid Rocket Engine Optimization for Small Launchers

Abstract

Hybrid rocket engines are a green alternative to solid rocket motors and may represent a low-cost alternative to kerosene fueled rockets, while granting performance and control features similar to that of typical storable liquid rocket engines. In this work, the design of a three-stage hybrid launcher is optimized by means of a coupled procedure: an evolutionary algorithm optimizes the engine design, whereas an indirect optimization method optimizes the corresponding ascent trajectory. The trajectory integration also provides the vertical emission profiles required for the evaluation of the environmental impact of the launch. The propellants are a paraffin-based wax and liquid oxygen. The vehicle is launched from the ground and uses an electric turbo pump feed system. The initial mass is given (5000 kg) and the insertion of the payload into a 600-km circular, and polar orbit is considered as a reference mission. Clusters of similar hybrid rocket engines, with only few differences, are employed in all stages to reduce the development and operational costs of the launcher. Optimization is carried out with the aim of maximizing the payload mass and then minimizing the overall environmental impact of the launch. The results show that satisfactory performance is achievable also considering rocket polluting emissions: the carbon footprint of the launch can be reduced by one fourth at the cost of a 5-kg payload mass reduction.

1. Introduction

In this century, the number of launches per year has grown steadily from about 50 successful launches in 2001 to more than 200 planned launches in 2022. The major drive of this positive trend is the shift of the launch industry from strictly governmental to a more commercial activity. In this scenario, rocket engine emissions in the atmosphere will increase at the same rate, gathering more and more attention from a regulatory point of view. To date, the amount of propellant burnt by rockets in the troposphere is about 0.01% of the fuel burnt annually by global aviation, and one can expect a proportional impact for their emissions. However, rockets are the main human-related source of localized pollutant emissions in the stratosphere and mesosphere, and hydrocarbon-based rocket engines are well known to produce far more carbon soot than conventional air-breathing engines [1].

There are two main architectures of rocket engines widely used in existing launchers: Liquid Rocket Engine (LRE) and Solid Rocket Motor (SRM). Their environmental impact depends primarily on the used propellants. SRMs usually consist of granules of a solid oxidizer, such as ammonium perchlorate and a polymeric binder with the addition of energetic additives, such as aluminum particles. Thus, their exhaust gases contain alumina particle and hydrochloric acid, which are well known to have a serious environmental impact.

LREs use several propellants combinations. Some common choices are represented by liquid oxygen (LOX)/kerosene (RP-1), liquid oxygen (LOX)/liquid hydrogen (LH2) and dinitrogen tetroxide (N₂O₄)/unsymmetrical dimethylhydrazine (UDMH). The exhaust gases of LOX/RP-1 engines contain water vapor, carbon dioxide and carbon soot. LOX/LH2 engines are employed more often in the upper stages, due to their high specific impulse, than in the first stage due to the low density of LH2. Unlike LOX/RP-1 engines, LOX/LH2 engines emit only water vapor. N₂O₄/UDMH are hypergolic storable propellants, which can guarantee high reliability, but their environmental impact is consistent due to the high toxicity of UDMH and the emission of nitric acid.

For this reason, toxic storables are seeing less and less application, while SRMs and kerosene-fueled LREs are still being widely used and will be employed also in the near future. Examples are the upcoming Ariane 6, which will be powered by the P120C solid rocket boosters [2] and SpaceX Merlin engines used on Falcon launchers [3].

Alongside the aforementioned heritage propulsion systems, the Hybrid Rocket Engine (HRE) is gathering more and more attention due to its good performance, low cost, and, which is the most crucial aspect in the present analysis, low environmental impact with respect to SRMs and storable LREs. These interesting benefits make hybrids’ competitors in some applications such as small launchers, as evidenced in recent years by the launch of Nammo’s Nucleus hybrid powered sounding rockets [4], ALTAIR, SMILE, and ENVOL European projects [5,6,7], and the three-stage small satellite launcher under development by HyImpulse [8]. However, in the past, some drawbacks hindered the spread and the exploitation of hybrid propulsion systems and slowed down their development, confining them still today to a lower technology readiness level with respect to LREs and SRMs. In particular, low regression rate and combustion instabilities can be regarded as the major issues [9]. Fortunately, in the last several decades, the research on liquefying fuels (e.g., paraffin-wax) and their unique entrainment phenomenon laid the foundation for the solution of the low regression rate problem [10], while the development of accurate numerical models allowed for the analysis of combustion instabilities [11]. For these reasons, in the present work, we consider paraffin-based wax and LOX as the propellants, which, in addition, do not emit alumina particles or other highly polluting compounds.

In the end, any combination of propellants causes some sort of emissions that may contain free radicals (e.g., NO, OH and Cl), inert (e.g., N₂ and CO₂) and under-oxidized compounds (e.g., H₂ and CO), resulting in indirect emissions formed downstream the rocket nozzles inside the high-temperature exhaust plume. Both direct and indirect emissions have been thoroughly analyzed for SRMs since the beginning of the Space Shuttle Program [12], due to concerns about stratospheric ozone depletion. On the contrary, to date, fewer analyses have been focused on the environmental impact of other rocket engines [13].

An additional critical point concerns the exact quantification of the environmental impact of a rocket engine, since there are many areas involved: stratospheric ozone depletion, greenhouse gases emissions in the troposphere, ecosystem toxicity, radiative forcing, mesospheric cloud formation, etc. A complete analysis is far beyond the scope of this paper, but the evaluation of all of the aforementioned environmental impacts must pass through the knowledge of direct emissions of a rocket engine during its operation.

In previous works, the authors analyzed the optimal design of hybrid rocket engines for several applications [14], including upper stages [15] and small satellite launchers [16]. In those works, the optimization aimed at maximizing the payload mass for a given target orbit, launcher configuration and initial mass. A coupled optimization approach was employed: a direct method optimized the engine design parameters, whereas the ascent trajectory was optimized by means of an indirect procedure. However, the emissions of the hybrid engines were not considered in the optimization merit function.

In this work, the authors embedded the numerical evaluation of the chemical composition of the exhausted gases inside the optimization procedure, assuming chemical equilibrium in the combustion chamber and frozen equilibrium expansion in the nozzle. In this way, the numerical procedure can precisely compute the vertical emission profile of the hybrid engines, and some environmental impact indices can be calculated (e.g., carbon footprint, instantaneous radiative forcing, etc.). These indices are used in the optimization procedure to minimize the environmental impact of the launch and obtain the corresponding optimized hybrid rocket engine design and ascent trajectory.

2. Hybrid Rocket Engine Modeling

The test case here considered is a hybrid-powered small satellite launcher. Its initial mass is equal to 5000 kg analogous to previous work by the authors [17]. The launcher consists of three stages, as other small launcher concepts, e.g., ALTAIR [5], SMILE [6] and SL1 by HyImpulse [18]. The number of almost identical HREs used in each stage (the few differences will be discussed in the following) is chosen to have a uniform split of the $\Delta v$, similar acceleration levels and mass fractions among the stages. Thus, a cluster of six engines is used in the first stage, three engines are used in the second stage, and a single HRE powers the third stage. LOX and a paraffin-based wax are the propellants, which grant good specific impulse and regression rate values. In addition, the high regression rate of the wax enables the use of a single port grain design and, at the same time, avoids the addition of aluminum particles in the fuel grain. The structural properties of the paraffin-wax grain may possibly be enhanced adding SEBS (styrene-ethylene-butylene-styrene block copolymer) in the fuel grain, without compromising the performance in terms of an obtained regression rate [19]. Hence, the proposed propellant combination does not emit alumina particles, unlike classical SRMs.

The HREs employ an electric turbo pump feed system, which consists of batteries, an electric motor and a pump. This solution has already been used in small-scale LREs, such as the Rutherford engine mounted on Rocket Lab’s Electron Launcher [20]. Moreover, the use of an electric pump is extremely appealing for hybrid application: oxidizer and fuel are stored in different physical states, and so it would be very difficult to generate the working fluid required to drive a turbine.

Commonality of engines drives the design, with the aim of reducing the overall cost of the launcher and the development time. Hence, some assumptions were made in order to use HREs as similar as possible in all the stages and at the same time match the constraint on exit pressure at liftoff, heat flux during ascent trajectory and maximum longitudinal acceleration at burnout. All fuel grains and oxidizer tanks are identical, but the third stage engine has a different oxidizer load because the feed system works at a lower power level to limit the longitudinal acceleration at the end of the ascent. The expansion area ratio of the nozzles is the same in the second and third stages, whereas, in the first stage, a smaller one is required to avoid separation issues at liftoff [21]. The mass of the interstage 1–2 adapter is assumed to be equal to 50 kg, the masses of the interstage 2–3 adapter and fairing (released at second stage jettisoning) sum up to 35 kg.

The engine modeling is based on the following assumptions and hypothesis:

1.

${\eta }_{c^{*}}=0.96$, correction factor to account for $c^{*}$ inefficiencies [22];

2.

${\eta }_F=0.98$, accounting for losses in the evaluation of the thrust coefficient $C_F$ [22];

3.

frozen equilibrium expansion, i.e., $\gamma =\mathrm{constant}$ throughout the nozzle;

4.

isentropic expansion in the nozzle [22];

5.

single circular port grain design, uniform regression rate along port axis, no residual fuel at engines burnout and no contribution of grain lateral end to overall combustion;

6.

regression rate correlation coefficient $a=9.1·{10}^{-5}{\mathrm{m}}^{2n+1}{\mathrm{s}}^{n-1}{\mathrm{kg}}^{-n}$ and exponent $n=0.69$ [10];

7.

third-order polynomial fittings of $c^{*}\left(\alpha \right)$ and $\gamma \left(\alpha \right)$ (computed with $p_c=10$ bar $=\mathrm{constant}$) [23];

8.

hydraulic resistance $Z=\mathrm{constant}$ and incompressible turbulent flow in the tank-chamber flow path;

9.

${\left(p_c\right)}_i=0.4{\left(p_d\right)}_i$ to avoid the coupling between the engine and the feed system during operation;

10.

$J=0.5$ at engines ignition to limit pressure losses in the combustion chamber;

11.

nozzle throat erosion is neglected;

12.

$ϵ=7$ for the first stage (to avoid separation at liftoff) and $ϵ=14$ for the upper stages, analogous to Ref. [24];

13.

fixed initial acceleration at liftoff (1.4 g), i.e., fixed initial thrust $F_i=11.5 \mathrm{kN}$ for each HRE, analogous to Ref. [24];

14.

oxidizer pump operation at constant power $P_E=\mathrm{constant}$;

15.

overall aspect ratio $L/D=10$ for the HREs, analogous to Ref. [24];

16.

ambient condition $p_0\left(h\right)$, ${\rho }_0\left(h\right)$ and $T_0\left(h\right)$ are known function of the altitude h;

17.

tank pressure $p_t=1 \mathrm{bar}=\mathrm{constant}$ in the oxidizer tank;

18.

negligible pressurizing gas mass required to have $p_t=\mathrm{constant}$ during engines operation.

Given this set of assumptions, two design parameters determine hybrid engine design and performance: the initial mixture ratio ${\alpha }_i$ and the initial pump discharge pressure ${\left(p_d\right)}_i$. Figure 1 reports a conceptual scheme of the considered HRE.

The specific heat ratio ${\gamma }_i\left({\alpha }_i\right)$ and characteristic velocity $c_i^{*}\left({\alpha }_i\right)$ at engine ignition are known. The initial chamber pressure is ${\left(p_c\right)}_i=0.4{\left(p_d\right)}_i$, and ${(p_e/p_c)}_i$ can be calculated from $ϵ$. The ambient pressure ${\left(p_0\right)}_i$ is known at launch location. Therefore, the thrust coefficient ${\left(C_F\right)}_i$ at liftoff can be calculated using Equation (1):

$$C_F={\eta }_F\left\{\sqrt{\frac{{\displaystyle 2{\gamma }^2}}{{\displaystyle \gamma -1}}{\left(\frac{{\displaystyle 2}}{{\displaystyle \gamma +1}}\right)}^{\frac{{\displaystyle \gamma +1}}{{\displaystyle \gamma -1}}}\left[1-{\left(\frac{{\displaystyle p_e}}{{\displaystyle p_c}}\right)}^{\frac{{\displaystyle \gamma -1}}{{\displaystyle \gamma }}}\right]}+ϵ\frac{p_e}{p_c}\right\}-ϵ\frac{p_0}{p_c}$$

(1)

The initial propellant mass flow rate, oxidizer mass flow rate and fuel flow rate are then:

$${\left({\dot{m}}_p\right)}_i=\frac{F_i}{c_i^{*}{\left(C_F\right)}_i}$$

(2)

$${\left({\dot{m}}_O\right)}_i=\frac{{\alpha }_i}{1+{\alpha }_i}{\left({\dot{m}}_p\right)}_i$$

(3)

$${\left({\dot{m}}_F\right)}_i=\frac{1}{1+{\alpha }_i}{\left({\dot{m}}_p\right)}_i$$

(4)

Equation (5) gives the throat area of the nozzle $A_{th}$:

$$A_{th}=\frac{{\left({\dot{m}}_p\right)}_i{\left(c^{*}\right)}_i}{{\left(p_c\right)}_i}$$

(5)

Obviously, the exit area of the nozzle is $A_e=ϵA_{th}$. Once $A_{th}$ is known, the initial port area ${\left(A_p\right)}_i$ and grain inner radius are:

$${\left(A_p\right)}_i=A_{th}/J$$

(6)

$$R_i=\sqrt{\frac{A_{th}}{J\pi }}$$

(7)

The initial fuel regression rate can be evaluated from ${\left(A_p\right)}_i$ and ${\left({\dot{m}}_O\right)}_i$ as:

$${\left(\dot{r}\right)}_i={\left[\frac{\mathrm{d}R}{\mathrm{d}t}\right]}_i=a{\left[\frac{{\left({\dot{m}}_O\right)}_i}{{\left(A_p\right)}_i}\right]}^n$$

(8)

The initial burning surface ${\left(A_b\right)}_i$ depends on ${\left(\dot{r}\right)}_i$, ${\rho }_F$ and ${\left({\dot{m}}_F\right)}_i$ as per:

$${\left(A_b\right)}_i=\frac{{\left({\dot{m}}_F\right)}_i}{{\rho }_F{\left(\dot{r}\right)}_i}$$

(9)

In the end, the grain length $l_g$ is:

$$l_g=\frac{{\left(A_b\right)}_i}{2\pi R_i}$$

(10)

At this point, the initial geometry of the grain and the nozzle, the engine performance and the mass flow rates are known parameters. Engine operation is ruled by the feed system characteristic, which in this case is provided by the following equation:

$$P_E=\frac{{\displaystyle {\dot{m}}_O(p_d-p_t)}}{{\displaystyle {\rho }_O{\eta }_{ep}}}=\mathrm{constant}$$

(11)

The constant power level $P_E$ is fixed by ${\left(p_d\right)}_i$ (which is a design parameter) and ${\left({\dot{m}}_O\right)}_i$ (known from Equation (3)).

The time integration of Equation (8) during operation gives the updated values of the grain radius R, the port area $A_p$ and the burning area $A_b$, which in turn gives the fuel flow rate ${\dot{m}}_F$ and the mixture ratio $\alpha$:

$$\alpha =\frac{{\dot{m}}_O}{{\dot{m}}_F}$$

(12)

Thus, during engine operation, $p_d$ is given by matching the pump and injector equations for ${\dot{m}}_O$, which is numerically solved by Newton’s method:

$${\dot{m}}_O=\frac{{\displaystyle P_E(p_d-p_t)}}{{\displaystyle {\rho }_O{\eta }_{ep}}}=\sqrt{\frac{p_d-p_1}{Z}}$$

(13)

where the overall efficiency ${\eta }_{ep}=0.53$ is introduced to account for the losses in the conversion process of electrical energy into flow head rise [25]. Chamber head end pressure $p_1$ depends on $p_c$ as per [26]:

$$p_1=\left[1+0.2{\left(\frac{A_{th}}{A_p}\right)}^2\right]p_c$$

(14)

and $p_c$ is:

$$p_c=\frac{({\dot{m}}_O+{\dot{m}}_F)c^{*}}{A_{th}}$$

(15)

where $c^{*}$ is computed for the actual mixture ratio $\alpha$ obtained from Equation (12).

The thrust level is $F=p_c{A}_{th}{C}_F$, where $C_F$ is evaluated at the actual altitude, in order to integrate the trajectory equations and for the actual value of $\gamma \left(\alpha \right)$. The integration of Equations (8)–(15) until engine burnout gives the grain outer radius $R_f$, the web thickness $w=R_f-R_i$, and the overall propellants mass exhausted $m_p=m_O+m_F$.

At this point, one can estimate the dry mass of the system $m_d$ as the overall mass of structural components (combustion chamber, nozzle, tank and rocket casing) plus electric components (motor, pump and battery). The structural masses are estimated as detailed in Ref. [16]. Specific power and energy density values (${\delta }_{ep}=3.92\mathrm{kW}/\mathrm{kg}$}, ${\delta }_{bp}=6.95\mathrm{kW}/\mathrm{kg}$} and ${\delta }_{be}=198.5\mathrm{Wh}/\mathrm{kg}$}) are introduced to compute the electric motor and pump mass $m_{ep}$ and the battery mass $m_b$, as follows [25]:

$$m_{ep}=\frac{{\displaystyle P_{E,max}}}{{\displaystyle {\delta }_{ep}}}$$

(16)

$$m_b=1.2max\left(\frac{{\displaystyle P_{E,max}}}{{\displaystyle {\delta }_{bp}}},\frac{{\displaystyle E_{E,tot}}}{{\displaystyle {\delta }_{be}}}\right)$$

(17)

The latter equation can be rewritten noting that $E_{E,tot}=P_E{t}_{burn}$ (because $P_E=\mathrm{constant}$) and introducing the characteristic burn time $t_{burn}^{*}={\delta }_{be}/{\delta }_{bp}=103\mathrm{s}$:

$$m_b=1.2\frac{{\displaystyle P_E}}{{\displaystyle {\delta }_{bp}}}max\left(1,\frac{{\displaystyle t_{burn}}}{{\displaystyle t_{burn}^{*}}}\right)$$

(18)

Hence, the power-constrained mass $m_b=1.2P_E/{\delta }_{bp}$ is known before trajectory evaluation (i.e., it only depends on the engine design) and represents the actual battery mass when $t_{burn}\le t_{burn}^{*}$ (i.e., short engine burns). On the contrary, if $t_{burn}\ge t_{burn}^{*}$ (i.e., long engine burns), the actual battery mass is given by the energy constraint and must be updated a posteriori (reducing the payload mass fraction) once the actual burn duration is known from trajectory integration.

The propulsion system dry mass is given by the sum of combustion chamber, nozzle, tank, rocket casing and electric components mass. The structural masses are estimated by means of suitable assumptions and approximations. A 6-mm insulating liner (with density equal to that of the solid fuel) and an aluminum alloy cylindrical wall encloses the combustion chamber. The aluminum oxidizer tank is cylindrical, with spherical ends, and its diameter is selected in order to obtain the chosen overall length-to-diameter ratio (here $L/D=10$). The wall thicknesses of tank and combustion chamber are selected to withstand internal pressure, assuming a 1.25 safety factor. In addition, a 1-mm thick cylindrical aluminum casing encapsulates each HRE. A 45-deg convergent and a 20-deg divergent nozzle (half-opening angle) with a phenolic silica ablative layer is taken into account. The nozzle mass is evaluated assuming a uniform thickness equal to half the value obtained according to Ref. [27] for the throat thickness; here, the authors employ average values of the transport properties and an estimation of the heat flux at ignition. In the end, nozzle structural mass is small with respect to the ablative layer mass and is neglected.

An in-house method performs the coupled optimization of engine design and ascent trajectory. A tentative engine design is assumed (i.e., ${\alpha }_i$ and ${\left(p_d\right)}_i$ are given), fixing the hybrid engine performance. Then, an indirect method optimizes the corresponding ascent trajectory for assigned target orbit and constraints (as will be discussed in Section 3) [28]. At engine burnout, propellant, dry and payload masses are computed, along with the environmental impact indices that will be described in Section 4.

At this point, the optimization merit function can be computed, and a particle swarm optimization (PSO) algorithm is used to optimize the engine design. Details about PSO algorithm are provided in Ref. [29];

Table 1. PSO settings.

Parameter	Value
Number of generations, $N_G$	100
Number of particles, $N_I$	20
Dimension of particles	2
Ranges of particles	${\mathit{b}}_U-{\mathit{b}}_L$
PSO method	1-trelea type 1
Cognitive acceleration, $C1$	2.0
Social acceleration, $C2$	2.0
Check population method	Saturation
End velocity weight	0.4
Linear varying factor	0.2
Maximum velocity, $v_{max}$	$0.25({\mathit{b}}_U-{\mathit{b}}_L)$
Mass mutation parameter	98%

and

Table 2. Design parameters’ ranges.

Boundary	${\mathit{\alpha }}_{\mathit{i}}$	${\left({\mathit{p}}_{\mathit{d}}\right)}_{\mathit{i}}$
	-	bar
${\mathit{b}}_L$	1.5	40
${\mathit{b}}_U$	2.5	60

show the setup and the parameters ranges used in the optimization. Figure 2 shows a simple flow-chart of the coupled optimization procedure.

3. Trajectory Optimization

In this work, a ground launch is considered from a 45-degree latitude. The target orbit is specified by means of altitude (600 km), inclination (95-deg) and eccentricity (zero, i.e., circular orbit). The authors chose this as a reference orbit analogous to the ENVOL project specification [7]. The trajectory is split in eight phases:

• Vertical ascent phase (fixed final altitude), which is followed by an instantaneous velocity rotation;
• Zero-lift gravity-turn ascent until first stage burnout (time duration is free and optimized);
• 8-second coast arc for first stage jettisoning;
• second stage burn (thrust direction and time duration are free and optimized);
• 8-second coast arc for second stage and payload fairing jettisoning;
• Third-stage burn (thrust direction and time duration are free and optimized);
• Coast arc (free and optimized length);
• Third-stage burn for orbit circularization (thrust direction and time duration are free and optimized).

The velocity rotation maneuver at the end of phase 1 usually has a duration of only a few seconds (e.g., 5), which is much shorter than the thrusting time and has a negligible impact on performance evaluation. Phase 8 tends to be short, and hence the circularization burn could eventually be performed by a payload engine, if available on board, instead of re-igniting third stage HRE. Re-ignition capability of liquefying fuels in space can be achieved using multiple ignition systems and adding SEBS (styrene-ethylene-butylene-styrene block copolymer) to the paraffin-based wax [19]. The free molecular heat flux $\varphi =(1/2){\rho }_{atm}{v}_{rel}^3$ has to remain below 1.136 $\mathrm{kW}/{\mathrm{m}}^2$} after the jettisoning of the payload fairing, which occurs during phase 5, to avoid thermal damage to the payload. Therefore, the ascent trajectory is forced to be steep in the early phases, imposing an altitude constraint at the first shutdown of the third stage (end of the sixth phase).

The trajectory optimization maximizes the payload mass delivered to the target orbit, given the engine design. A point-mass rocket is considered and the state equations in an inertial Earth-centered reference frame are:

$$\frac{\mathrm{d}\mathit{r}}{\mathrm{d}t}=\mathit{v}\frac{\mathrm{d}\mathit{v}}{\mathrm{d}t}=-G{M}_{\oplus }\frac{\left|\right|{\mathit{r}}^3\left|\right|}{\mathit{r}}+\frac{\mathit{F}-\mathit{D}}{m}\frac{\mathrm{d}m}{\mathrm{d}t}=-\frac{F}{c}$$

(19)

The authors assumed an inverse square gravity field and wrote the state equation in non-dimensional form to obtain an high precision in their numerical integration. The aerodynamic drag is $D=(1/2){\rho }_{atm}{C}_D{S}_{ref}{v}_{rel}^2$, where $S_{ref}$ is the sum of the cross sections of the HREs used in each stage, and ${\mathit{v}}_{rel}=\mathit{v}-\mathbf{\omega }\times \mathit{r}$ ($\mathbf{\omega }$ is the Earth’s angular velocity). The drag coefficient $C_D$ depends on the Mach number as per the typical rockets law shown in Figure 3.

In the indirect trajectory optimization procedure, the vacuum thrust $F_{vac}$ is a known function of time because the engine design is given by the design parameters. Therefore, the thrust F is:

$$F=F_{vac}-ϵA_{th}{p}_{atm}$$

(20)

where ambient condition $p_{amb}\left(h\right)$ (and also $T_{amb}\left(h\right)$ and ${\rho }_{amb}\left(h\right)$) are known as functions of altitude h in accordance to numerical fits of Earth’s atmosphere.

The boundary conditions at the initial time $t_0=0$ fix the initial values of the state variables (i.e., position, velocity and mass). Altitude, inclination and eccentricity of the target orbit give the boundary conditions at the final time $t_f$, fixing the final radius, vertical velocity, horizontal velocity magnitude and inclination. Although the HREs employed in each stage are similar, they are not identical. The constraints on burning times (i.e., $t_2=t_4-t_3=t_8-t_7+t_6-t_5$) must be replaced by mass constraints on the exhausted propellant masses. Hence, aiming at the highest commonality, the constraints are imposed on the exhausted fuel masses, i.e., all the fuel grains are identical. In the third stage, less oxidizer is required (due to the different power level) and the corresponding tank can be filled only partially, reserving the saved mass for the payload. The reformulation of these boundary conditions for optimality is here omitted for the sake of conciseness. The duration of the coast arcs for the stage jettisoning is given and $t_3-t_2=t_5-t_4$. The altitude at $t_6$ (end of the third burn) is fixed to avoid excessive heat flux after fairing jettisoning.

The Optimal Control Theory (OCT) optimizes the ascent trajectory, once the performance of each stage is given. First, adjoint variables $\lambda$ are associated with the problem state equations. Then, the Hamiltonian, whose formal expression depends on the phase of flight, is defined as reported in Equation (21):

$$H={\mathit{\lambda }}_r\mathit{v}+{\mathit{\lambda }}_v\left(\frac{\mathit{r}}{{\left|\mathit{r}\right|}^3}+\frac{\mathit{F}-\mathit{D}}{m}\right)-{\lambda }_m\frac{F}{c}$$

(21)

The OCT provides the Euler–Lagrange equations for the adjoint variables:

$$\frac{\mathrm{d}{\mathit{\lambda }}_r}{\mathrm{d}t}=-\frac{\mathrm{d}H}{\mathrm{d}\mathit{r}}\frac{\mathrm{d}{\mathit{\lambda }}_v}{\mathrm{d}t}=-\frac{\mathrm{d}H}{\mathrm{d}\mathit{v}}\frac{\mathrm{d}{\lambda }_m}{\mathrm{d}t}=-\frac{\mathrm{d}H}{\mathrm{d}m}$$

(22)

In the considered mission, the thrust is vertical during phase (1), parallel to the relative velocity vector during phase (2), free and optimized throughout phases (4), (6) and (8), and null during phases (3), (5) and (7) (coast arcs). The optimal thrust direction during phases (4), (6) and (8) results in being parallel to the velocity adjoint vector, also known as the primer vector. The OCT also provides the boundary conditions for optimality at the initial and final points, and at the boundary points of each phase of flight.

The derivation of the optimality conditions is easier for the dual problem of maximization of the insertion radius for a given HRE mass (i.e., given payload). Velocity adjoint vector and relative velocity vector must be parallel at $t_1$. At the final time $t_f$, position and velocity values are related to their adjoint variables by two additional conditions (here omitted). The OCT provides the transversality conditions to determine the relevant times. Time is formally free in this formulation because time boundaries become mass constraints. Hence, the Hamiltonian must be continuous at $t_6$ and $t_7$ and null at $t_8$. The radius adjoint variable has a free discontinuity at $t_6$, which is an additional optimization variable. The condition ${\lambda }_{rf}=1$ is replaced by ${\lambda }_{r0}=1$. The use of the OCT leads to a multipoint boundary value problem (BVP), which is solved by an iterative numerical procedure based on Newton’s method [28].

4. Emissions Analysis

Rocket emission refers to the vertical distribution and composition of the exhausted gases, which obviously depends on the propellants choice. The fuel consumption of aviation out-scales the rocket propellant consumption in the troposphere by a factor of about ${10}^2$, but rockets represent the main source of direct human-related emissions above 15 km altitude, i.e., in the stratosphere and in the mesosphere. The lifetime of combustion products in the troposphere is limited and and their presence can be neglected, with the exception of carbon dioxide. In the stratosphere, the lifetime of combustion products is longer, and the distribution of four of them is usually analyzed: carbon dioxide (CO₂), water vapor (H₂O), black carbon (BC), and submicrometer alumina particles (SMF) [1].

The conventional approach to emissions quantification is based on some simplifications. Emission Indices (EIs) are synthetic factors, expressed in grams per kg of burned propellant, which can be used to estimate the steady state burden M for a specific compound as per the following equation:

$$M=\frac{2}{3}PN{10}^{3}EI\tau$$

(23)

where the factor $2/3$ accounts for typical ground launches (that is, two-thirds of the propellants are exhausted above tropopause), P is the total propellant mass of the launcher expressed in tons, N is the number of launches per year, $\tau =4$ years is the average stratospheric lifetime, and the 10${}^3$ factor is introduced to obtain the steady state burden M in kilograms. Emission indices do not take afterburning in the rocket plume into account, since it usually occurs at low altitude. Typical EIs for hybrids from the literature are: EI(CO${}_2$) = 200, EI(H₂O) = 200, EI(BC) = 40, EI(SMF) = 0 [1].

In order to overcome these assumptions, in this work, the authors embed NASA CEA (Chemical Equilibrium with Application) [23] code in the trajectory optimization to precisely compute the vertical emission profile of the three stage launcher here considered. In this way, no assumptions are needed on the emissions distribution throughout the ascent, and the use of EIs is avoided. The simulations run in “rocket” mode and the input parameters are chamber pressure $p_c$ and mixture ratio $\alpha$ calculated by the engine model throughout the trajectory integration. The oxidizer is liquid oxygen, while the paraffin-based fuel is specified by its gross formula (C₃₂H₆₆) and enthalpy of formation ($\Delta H_f=-697.2\mathrm{kJ}/\mathrm{mol}$} at 533 K) [30]. Frozen equilibrium expansion in the nozzle is assumed for the simulations.

Preliminary results suggested that carbon soot formation occurs, for the considered propellants, when $\alpha \le 1.15$ for chamber pressures in the range 20 bar $\le p_c\le$ 60 bar (see Figure 4), following the C/O=1 rule for the thermodynamic onset of soot production. Therefore, despite the lack of a kinetic model, carbon soot formation can be neglected in the proposed approach, since the optimal mixture ratio (that is, the mixture ratio that maximizes the specific impulse) is around 2 for the considered propellant combination [31].

Therefore, the following combustion products are taken into account: carbon monoxide (CO), carbon dioxide (CO${}_2$), monoatomic hydrogen (H), gaseous hydrogen (H${}_2$), water (H${}_2$O), monoatomic oxygen (O), gaseous oxygen (O${}_2$), hydroxide (OH), ethylene (C${}_2$H₄) and methane (CH₄) [32].

This numerical procedure provides the mass fractions of the combustion products as a function of the rocket altitude $X_i\left(h\left(t\right)\right)$, which can be integrated between 17 km (average tropopause altitude) and 50 km (average stratopause altitude) to evaluate the actual stratospheric EIs as per:

$$EI\left(C{O}_2\right)={\int }_{t_T}^{t_S}{X}_{C{O}_2}\left(h\left(t\right)\right){\dot{m}}_p\left(t\right)dt$$

(24)

$$EI\left(H_{2}O\right)={\int }_{t_T}^{t_S}{X}_{H_{2}O}\left(h\left(t\right)\right){\dot{m}}_p\left(t\right)dt$$

(25)

with $h\left(t_T\right)=17\mathrm{km}$ and $h\left(t_S\right)=50\mathrm{km}$. The emission index of black carbon EI(BC) = 0 because carbon soot formation is excluded in this analysis as discussed above, and the emission index of alumina particles EI(SMF) = 0 too because the paraffin-based wax avoids the presence of aluminum in the fuel grain. The knowledge of the exact vertical emission profile allows also a more precise evaluation of the steady state burden M for a specific component. In fact, the term $\frac{2}{3}P{10}^{3}EI$, which approximates the launcher emissions in the troposphere as a function of the overall propellant mass P, can be easily replaced by the actual mass of the exhausted component obtained by integration. In particular, the considered hybrid launcher does not follow the “two third rule”. In fact, in preliminary optimization runs, around 25% of the propellant burns in the troposphere and around 75% burns above the tropopause, following a “three fourth rule”.

The environmental impact of rockets can be analyzed from a different point of view by evaluating the instantaneous Radiative Forcing (iRF) from specific exhausted gases. The iRF evaluates the immediate change in the heat flux at the top of the atmosphere after the insertion of each chemical component into a model stratosphere. The approach to evaluate the iRF is different for each component, as detailed in Ref. [1].

The iRF from CO${}_2$ can be calculated by means of Equation (26):

$$iR{F}_{C{O}_2}=l_{LW}{\sigma }_{C{O}_2}N{M}_{C{O}_2}{A}_E^{-1}$$

(26)

where $l_{LW}=235\mathrm{W}/{\mathrm{m}}^2$ is the terrestrial long wave integrated flux (wave length from $6.25·{10}^{-6}\mathrm{m}$ to $30·{10}^{-6}\mathrm{m}$), ${\sigma }_{C{O}_2}=0.005 {\mathrm{m}}^2/\mathrm{kg}$ is the mass-specific absorption coefficient of CO${}_2$, $N=1$ to evaluate the iRF from a single launch, $M_{C{O}_2}$ is the steady state stratospheric burden for CO${}_2$, and $A_E$ is the solar illuminated area of the surface of the Earth.

The emission of water vapor in the stratosphere warms the troposphere (long-wave and short-wave absorption) while cooling the stratosphere (long-wave emission). The iRF from H${}_2$O is calculated by means of Equation (27):

$$iR{F}_{H_{2}O}=l_{LW}{\sigma }_{H_{2}O}N{M}_{H_{2}O}{A}^{-1}$$

(27)

where ${\sigma }_{H_{2}O}=4 {\mathrm{m}}^2/\mathrm{kg}$ is the mass-specific absorption coefficient of H${}_2$O, $M_{C{O}_2}$ is the steady state stratospheric burden for H${}_2$O, and $A=1.2·{10}^{14} {\mathrm{m}}^2$ is the surface area of the accumulation region for exhausted gases, which extends between 15 and 30 km altitude and 30-deg N-80-deg N latitude, as suggested by global climate models [33].

In the general case, the radiative forcing from particles, both BC and alumina, is of concern for rockets because their accumulation in the stratosphere scatters and absorbs the solar short wave radiation, while absorbing terrestrial long wave up-welling radiation. These phenomena cause a localized heating of the stratosphere, which contributes to ozone depletion. However, the considered HREs do not produce alumina or BC particles, as discussed above, hence their iRF is not analyzed.

In the end, the knowledge of the composition of exhaust gases in the troposphere can be used to evaluate the Carbon FootPrint (CFP) of a single launch. Global Warming Potentials (GWP) relative to CO${}_2$ are defined for the most common greenhouse gases, and express the amount of heat absorbed by a specific greenhouse compound in the atmosphere, as a multiple of the heat absorbed by the same mass of CO${}_2$ (i.e., $GW{P}_{C{O}_2}=1$) [34,35].

Table 3. GWP (lifetime 100 years).

Gas	GWP
CO${}_2$	1
CO	3
CH${}_4$	27.9
H${}_2$	5.8
H	5.8

reports the GWP for the gases of interested for this analysis.

The carbon footprint of a single launch can be calculated as:

$$CFP=\sum _i^n\left(GW{P}_i{\int }_{t_0}^{t_f}{\dot{m}}_p\left(t\right)X_i\left(h\left(t\right)\right)dt\right)$$

(28)

where $GW{P}_i$ and $X_i\left(h\left(t\right)\right)$ are the global warming potential and the mass fraction of the i-th combustion product. In this analysis, five greenhouse gases are taken into account (that is, $n=5$): carbon monoxide (CO), carbon dioxide (CO${}_2$), methane (CH${}_4$), gaseous hydrogen (H${}_2$) and monoatomic hydrogen (H). CFP does not take into account water (H${}_2$O), monoatomic oxygen (O), gaseous oxygen (O${}_2$) and hydroxide (OH) because their GWP is zero; data about ethylene can not be found in the literature so its GWP is not considered here. In any case, the hypothesis of chemical equilibrium leads to very small mass fractions for both methane and ethylene during the ascent, with respect to CO₂ and CO fractions, making negligible their possible contribution to CFP.

The following Section 5 reports vertical emission profiles, EIs, iRFs and CFPs of optimized designs of the proposed hybrid launcher.

5. Results

First, the design of the proposed launcher is optimized aiming at the maximization of the payload delivered in the specified orbit (600-km circular, polar orbit). The preliminary results showed that an unconstrained optimization leads to unfeasible solutions. The launcher payload mass benefits from values of ${\left(p_d\right)}_i$ in the upper end of the design range (i.e., ${\left(p_d\right)}_i\ge 50 \mathrm{bar}$) because of the reduction of the engine dry mass: the initial thrust $F_i$ of the HREs is fixed; thus, higher ${\left(p_c\right)}_i=0.4{\left(p_d\right)}_i$ results in smaller $A_t$, which in turn makes the nozzle lighter since $ϵ$ is given too. At the same time, the increase of ${\left(p_d\right)}_i$ makes the electric motor and pump heavier (as per Equation (11) and 16), which partially compensates for the overall reduction in dry mass of the engine. This trend is beneficial in terms of payload mass, but unfortunately leads to excessive longitudinal acceleration at stages’ burnouts. For this reason, a constraint is imposed in the optimization procedure that forces the maximum longitudinal acceleration experienced by the launcher below 6.5 g. Thus, the optimization problem can be formally cast as:

$$\begin{array}{cccc}& \mathrm{find}\hfill & & \mathit{b}\in {\mathbb{R}}^n\hfill \\ & \mathrm{to}\mathrm{maximize}\hfill & & \Phi \left(\mathit{b}\right)\hfill \\ & \mathrm{subject}\mathrm{to}\hfill & & g_j\left(\mathit{b}\right)⩽0,j=1,\dots ,r\hfill \\ & \mathrm{and}\mathrm{to}\hfill & & {\mathit{b}}_L⩽\mathit{b}⩽{\mathit{b}}_U\hfill \end{array}$$

(29)

where $\Phi \left(\mathit{b}\right)$ is here the launcher payload mass $\mu$, $\mathit{b}=[{\alpha }_i,{\left(p_d\right)}_i]$ is the optimization parameters vector, $g_j$ is the j-th inequality constraint imposed in the trajectory optimization (see Section 3), and $b_L$ and $b_U$ are the lower and upper boundary of the design parameters as reported in Table 2.

Table 4. Comparison of environmental impact indices. EIs units are grams emitted per kilogram of propellant burned. $CFP/\mu$ units are kg of CO${}_2$ equivalent emitted per kg of payload mass delivered in the target orbit.

$\mathit{C}\mathit{a}\mathit{s}\mathit{e}$	${\mathit{\alpha }}_{\mathit{i}}$	${\left({\mathit{p}}_{\mathit{d}}\right)}_{\mathit{i}}$	$\mathit{\mu }$	$\mathit{C}\mathit{F}\mathit{P}$	$\frac{\mathit{C}\mathit{F}\mathit{P}}{\mathit{\mu }}$	$\mathit{E}\mathit{I}\left(\mathit{C}{\mathit{O}}_2\right)$	$\mathit{E}\mathit{I}\left({\mathit{H}}_2\mathit{O}\right)$	$\mathit{i}\mathit{R}{\mathit{F}}_{\mathit{C}{\mathit{O}}_2}$	$\mathit{i}\mathit{R}{\mathit{F}}_{{\mathit{H}}_2\mathit{O}}$
	-	bar	kg	kg	-	-	-	mW/m${}^2$	mW/m${}^2$
REF	1.81	50.69	51.1	6636	130.0	278.9	254.5	4.944$·{10}^{-12}$	1.023$·{10}^{-8}$
OPT	2.33	50.48	44.8	4968	111.0	370.1	251.1	5.973$·{10}^{-12}$	9.185$·{10}^{-9}$

shows the optimized design and performance as REF, alongside the corresponding environmental impact indices evaluated as described in Section 4. For the sake of completeness, the following Figure 5 reports $p_c$, $\alpha$ and ${\dot{m}}_p$ vs. launcher altitude h during the ascent up to 100 km. In this way, interested readers can perform different CEA simulations or emissions evaluations using other models for any given launcher altitude below 100 km.

Figure 6 shows the relation between the mass fractions of the considered exhausted gases and the rocket altitude in the lower atmosphere ($h\le 100\mathrm{km}$) for the REF optimized design. As expected, the mass fractions of C${}_2$H${}_4$ and CH${}_4$ are negligible (omitted in (b)). Since $\alpha$ and $p_c$ vary during HRE burns, the mass fractions of the exhausted gases are not constant. In particular, one can notice that among the major fractions (on the left), $X_{C{O}_2}$ and $X_{H_{2}O}$ increase, while $X_{CO}$ decreases due to the increase of $\alpha$ during the flight, which is a well-known phenomenon typical of hybrids called mixture ratio shifting. This behavior may be useful for managing, and we hope to reduce the overall emissions of engines in terms of CFP since $GW{P}_{CO}=3$, while $GP{W}_{C{O}_2}=1$ and $GW{P}_{H_{2}O}=0$. Analogously, Figure 7 shows the integral of rocket emission for each exhausted gas in the lower atmosphere ($h\le 100\mathrm{km}$), computed as:

$$m_i={\int }_{t_0}^{t_f}{\dot{m}}_p\left(t\right)X_i\left(h\left(t\right)\right)dt$$

(30)

These vertical emission profiles and indices establish a baseline level of environmental impact for the proposed hybrid-powered launcher. In order to minimize this impact, the optimization procedure is repeated switching to an emission-oriented merit function in place of the payload mass delivered. Possible choices are the CFP, the total iRF (i.e., $iR{F}_{C{O}_2}+iR{F}_{H_{2}O}$), and the payload mass specific CFP (i.e., $CFP/\mu$). However, preliminary results showed that the optimal designs obtained are really close to each other regardless of the actual emission index used in the optimization, all of them being functions of the overall rocket emission. Thus, in this work, the author selected the payload mass-specific CFP as a merit function to be minimized by the optimization procedure, and the optimization problem can be formally cast as:

$$\begin{array}{cccc}& \mathrm{find}\hfill & & \mathit{b}\in {\mathbb{R}}^n\hfill \\ & \mathrm{to}\mathrm{minimize}\hfill & & \Phi \left(\mathit{b}\right)\hfill \\ & \mathrm{subject}\mathrm{to}\hfill & & g_j\left(\mathit{b}\right)⩽0,j=1,\dots ,r\hfill \\ & \mathrm{and}\mathrm{to}\hfill & & {\mathit{b}}_L⩽\mathit{b}⩽{\mathit{b}}_U\hfill \end{array}$$

(31)

where $\Phi \left(\mathit{b}\right)$ is here $CFP/\mu$, $\mathit{b}=[{\alpha }_i,{\left(p_d\right)}_i]$ is the vector of optimization parameters, $g_j$ is the j-th inequality constraint imposed in the trajectory optimization (see Section 3) and ${\mathit{b}}_L$ and ${\mathit{b}}_U$ are the lower and upper boundary of the design parameters as reported in Table 2. The emission-driven optimized design is reported in Table 4 as OPT, alongside the corresponding environmental impact indices.

Analogously to the REF case, the value of ${\left(p_d\right)}_i$ is bounded by the longitudinal acceleration constraint for a given value of ${\alpha }_i$. The optimization procedure increases the value of ${\alpha }_i$ (from 1.81 to 2.33) in order to obtain the aforementioned reduction of $X_{CO}$ (red curves in Figure 8), which results in a relevant reduction of the CFP for the OPT solution (approximately 3/4 of the REF case), despite the increase of $X_{C{O}_2}$ (blue curves in Figure 8). At the same time, this choice of ${\alpha }_i$ is not optimal in terms of payload and the OPT solution obviously has payload mass smaller by roughly 12% with respect to the REF optimized solution (44.8 kg vs. 51.1 kg). However, the % reduction of CFP is twice the payload mass reduction, making the minimum $CFP/\mu$ design and the minimum $CFP$ design actually the same in this application. Looking at the $iRFs$ values reported in the Table 4, one can observe that $iR{F}_{H_{2}O}$ is three orders of magnitude higher than $iR{F}_{C{O}_2}$ for both the REF and the OPT solution. Thus, the contribution of $C{O}_2$ emissions to the total instantaneous radiative forcing can be neglected for the hybrid launcher analyzed. The increase of ${\alpha }_i$ of the OPT solution actually reduces the overall $H_{2}O$ emission in the stratosphere (although increases $C{O}_2$ emission) with respect to REF solution, as shown in Figure 8 (on the right, green curves), resulting in a smaller $iR{F}_{H_{2}O}$. For this reason, the design which minimizes the CFP is also the same that minimizes the total $iRF$ for the considered application.

In the end, looking at the EIs reported in the Table 4, the reader can see that all the computed values of EIs are quite different from the value provided in the literature for hybrid, which are EI(CO${}_2$) = EI(H${}_2$O) = 200 [1]. This mismatch may be due to the peculiar shifting of the mixture ratio of hybrids during operation, which causes a variation of the relative contribution of each exhausted gases during the ascent of the launcher, depending on the choice of the initial mixture ratio ${\alpha }_i$. Moreover, this analysis suggests that a comparison of the environmental impact of different hybrid engine designs cannot be achieved by means of EIs, since the computed value of EI(CO${}_2$) is actually greater for the OPT solution, with respect to the REF solution, despite the smaller value obtained for CFP and CFP/$\mu$.

6. Conclusions

The evaluation of possible emissions reduction and its impact on the launcher payload has been carried out using the solution that optimizes the payload (delivered to a reference target 600 km circular and polar orbit) as a reference. This optimal design has a payload mass slightly greater than 50 kg in accordance with previous works on the same topic, and the trajectory integration provides the vertical emission profiles during the ascent and the environmental impact of the launch by means of specific indices (carbon footprint and instantaneous radiative forcing).

The optimization procedure is then repeated with the aim of minimizing the environmental impact of the launch. The optimal design results in being the same regardless of the choice of the impact index to be minimized, granting the reduction of both the carbon footprint of the launch and the total instantaneous radiative forcing. Two important results are found from this analysis. First, emissions can be significantly reduced (−25% carbon footprint, −10% total instantaneous radiative forcing), at the cost of a 5-kg payload mass penalty (−10%). Second, the results show that the conventional approach to emission quantification, based on emission indices, is not well suited to describe hybrid emissions.