# Integrated Aerothermal-Flight Optimization of a Solid-fuel Ducted Rocket Scramjet Vehicle

Megha Bhardwaj, Subramanian Rajashankar, Narayan Ananthkrishnan,
Anurag Sharma, Wan-Do Ki, Hyuck-Joon Namkoung

► To cite this version:

Megha Bhardwaj, Subramanian Rajashankar, Narayan Ananthkrishnan, Anurag Sharma, Wan-Do Ki, et al.. Integrated Aerothermal-Flight Optimization of a Solid-fuel Ducted Rocket Scramjet Vehicle. 2025. [hal-04946366](https://hal.science/hal-04946366v1)

HAL Id: hal-04946366

https://hal.science/hal-04946366v1

Preprint submitted on 13 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Integrated Aerothermal-Flight Optimization of a Solid-fuel
Ducted Rocket Scramjet Vehicle

Megha Bhardwaj*, Subramanian Rajashankar†, Narayan Ananthkrishnan‡, and Anurag Sharma§
Yanxiki Tech, Pune, Maharashtra, 411001, India

Wan-Do Ki¶, Hyuck-Joon Namkoung||
Hyundai Rotem, Uiwang-si, Gyeonggi-do, Korea

The problem of optimizing the operation of a solid-fuel ducted rocket scramjet (DRSJ) vehicle to maximize its mission range has been taken up in this work. An integrated aerothermal-flight model, including the scramjet engine itself, the gas generator, and the flight trajectory, has been constructed for this purpose. The climb angle and the gas generator pressure at the start of climb are the chosen design parameters. Their variation throughout the mission is calculated by using a dynamic inversion scheme. There is a narrow band of values of the design parameters where the requirements of adequate thrust for climb as well as near-steady state cruise are met. The optimal solution is seen to lie at the boundary defined by two constraints —thermal choking in the scram combustor and limiting air mass flow rate into the engine intake. Time simulation of the mission profile with the selected values of the design parameters is conducted to verify the results.

## Nomenclature

$A$ = Duct area, m²

$A_b$ = Propellant burn area, m²

$A_C$ = Isolator core flow area, m²

$A^*, A_t$ = Gas generator nozzle throat area, m²

$A_w$ = Wetted area, m²

$C_L, C_D$ = Lift, drag coefficient, respectively

$C_f$ = Skin friction coefficient

$c_r, c_t$ = Root, tip chord, respectively, m

<footnote>*R & D Engineer, Yanxiki Tech, megha@yanxiki.com.</footnote>

<footnote>†Senior R & D Engineer, Yanxiki Tech, rajashankar@yanxiki.com.</footnote>

<footnote>‡Technical Consultant, Yanxiki Tech, akn.korea.19@gmail.com, akn@aero.iitb.ac.in. Associate Fellow, AIAA.</footnote>

<footnote>§Principal Engineer, Yanxiki Tech, anurag@yanxiki.com.</footnote>

<footnote>¶Senior Research Engineer, Guided Munitions Team, Hyundai Rotem, Korea, gongdol@hyundai-rotem.co.kr.</footnote>

<footnote>||Team Leader, Guided Munitions Team, Hyundai Rotem, Korea, namkoung@hyundai-rotem.co.kr.</footnote>
$D$ = Combustor exit height, m

$D_H$ = Hydraulic diameter, m

$dx$ = Distance along isolator duct

$f$ = fuel-to-air ratio

$g$ = Acceleration due to gravity, m/s²

$h$ = Altitude, m

$h_c, h_p, h_r$ = Enthalpy of combustion, products, reactants, respectively, J/kg

$M$ = Mach number

$\dot{m}$ = Mass flow rate, kg/s

$P$ = Pressure, Pa

$\dot{Q}_c$ = Rate of heat release, J/s

$Q, Q_{rel}$ = Heat release per unit mass, J/kg

$\dot{Q}_w$ = Rate of heat loss, J/s

$\dot{r}_{burn}$ = Propellant burn rate, m/s

$S_{ref}$ = Reference area, m²

$s$ = Wing span (reference length), m

$T$ = Thrust, N

$T_{ax}$ = Axial component of thrust, N

$T$ = Temperature, K

$T_t, T_0$ = Total temperature, K

$U$ = Flow velocity, m/s

$V$ = Vehicle inertial velocity, m/s

$V_c$ = Gas generator free volume, m³

$W$ = Vehicle weight, kg

$\alpha$ = Angle of attack, deg

$\beta$ = Shock angle, deg

$\gamma$ = Climb angle, deg

$\gamma_{air}$ = Ratio of specific heats of air

$\gamma_c$ = Ratio of specific heats of pyrolysis products

$\zeta_d, \zeta_F, \zeta_v$ = Correction factor for discharge, thrust, velocity, respectively

$\theta$ = Flow turn angle, deg

$v(M)$ = Prandtl-Meyer function
$\rho$ = Density, kg/m³

$\omega, \omega_p$ = Bandwidth, /s

=

Subscripts =

$a$ = air

$act$ = actual value

$b$ = burn, solid propellant

$c$ = gas generator chamber

$e$ = exit

$est$ = estimated value

$f$ = fuel

$g$ = gas

$i$ = injected

$id$ = ideal value

$s$ = solid

SET = set value

# I. Introduction

AIR-LAUNCHED or shore-based high-speed anti-ship missiles can be a potent threat to aircraft carriers and other naval battle groups. Present-generation supersonic missiles, like the Brahmos <fnref n="1" />, fueled by a liquid ramjet engine, can fly at speeds in excess of Mach 3 and attack targets at a range of several hundred kilometers. Their next-generation counterparts are expected to fly faster at hypersonic speeds and achieve a greater range of the order of a thousand kilometers. These missiles are likely to be powered by solid-fuel scramjet engines. It is therefore of interest to consider whether solid-fuel scramjet engines can be conceptualized to meet the mission requirements of such an engagement.

A solid-fuel scramjet engine is conceptually no different from a liquid-fueled one, except for the manner in which the solid fuel is stored and combusted. In general, three possible arrangements of the solid fuel are possible for use with a scramjet engine, as illustrated in Fig. 1:

1. The plain solid-fuel scramjet (SFSJ), sketched in Fig. 1a, represents the simplest and earliest version of solid fuel use in a scramjet mode. The solid fuel is placed along the walls of the scram combustor, making it a compact arrangement. The incoming air pyrolizes the solid fuel, which further combusts in the scram combustor. Obviously, regulating the fuel burn rate, and hence the thrust, is a problem. Besides, the fuel may not burn

<footnote n="1">¹https://www.brahmos.com/</footnote>
<figure><img src="image_2.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(c)</figcaption></figure>

Fig. 1 Schematic of different types of solid-fuel scramjets: (a) Solid-fuel scramjet (SFSJ) [1], (b) Solid-fuel dual-combustor scramjet (DCSJ) [2], (c) Solid-fuel ducted rocket scramjet (DRSJ)[3].

completely within the length of the combustor, making this arrangement somewhat inefficient; hence, this type is not very popular at present.

2. Solid-fuel dual-combustor scramjet (DCSJ) is another possible scheme, as shown in Fig. 1b. In this case, there are two combustors mounted in a tandem arrangement (hence the name 'dual combustor'). The forward (or primary) one is a subsonic ram combustor, which functions as a gas generator. The solid fuel is again placed along the walls of this ram combustor and is partially combusted with the air taken in from a separate ram intake. Then, the fuel-rich gas from the primary combustor is injected into the secondary scram combustor, where it mixes with air from a scram intake and further combustion takes place. While this arrangement is favorable for complete combustion of the solid fuel, regulating the fuel burn rate, and hence the thrust, remains an issue. The need for distinct intakes makes the intake design a tricky affair, with the ram intake, which slows down the hypersonic free stream flow to a low subsonic speed, suffering a very high loss of total pressure. Some control is possible in
<figure><img src="image_5.png" /><figcaption>Fig. 2 Typical waverider configuration for a ducted rocket scramjet vehicle [11].</figcaption></figure>

the rate of injection of fuel-rich gas from the exit of the ram combustor. However, a build-up of pressure in the ram combustor may result in the normal shock system in the ram intake being ejected forward out of the intake duct causing *unstart* from which recovery may be difficult. The axisymmetric geometry is often not favored for hypersonic cruise flight. The tandem combustor arrangement tends to increase the vehicle length, though the scram combustor may possibly be made shorter in this case thanks to the use of hot, gaseous fuel being injected. These significant design issues have made this arrangement less popular.

3. Solid-fuel ducted rocket scramjet (DRSJ), as indicated in Fig. 1c, consists of a gas generator loaded with a fuel-rich solid propellant. On combustion, a hot, fuel-rich gaseous fuel is created, which is introduced into the scram combustor, where it burns relatively easily. The gas generator in the DRSJ does not require atmospheric oxygen, hence a single intake for the scram combustor suffices for this configuration. By controlling the outflow from the gas generator to the scram combustor, a degree of throttling is possible allowing for control of thrust. Thus, the DRSJ overcomes many of the shortcomings of the other two solid-fuel scramjet schemes discussed above.

Of these, the ducted rocket scramjet (DRSJ) is easily the most promising configuration for a solid-fuel-powered hypersonic cruise vehicle. Its concept was initially proposed by Lv et al. [3, 4], and further studies demonstrated that it was viable for the DRSJ to operate at a high altitude with a low ambient pressure [1]. Following that, there has been heightened research interest towards this configuration [5–9]. The use of metal additives with the solid propellant in DRSJ gas generators to boost the thrust has been studied [10]. Particles of Boron, Magnesium, Aluminum embedded in the solid propellant grain appear to provide higher equilibrium combustion temperatures compared to pure hydrocarbon-based fuels. However, these particles do not usually burn completely in the gas generator and are injected into the scram combustor in solid phase along with the gaseous fuel where they need to combust further.

For a hypersonic cruise mission, it is common to consider a DRSJ vehicle with a waverider configuration, as illustrated in Fig. 2. In this case, the engine flow path is integral with the vehicle under-body. The waverider is
<figure><img src="image_6.png" /><figcaption>**Fig. 3 Three-phase mission profile for ducted rocket scramjet considered in this study.**</figcaption></figure>

characterized by high lift-to-drag ratio (L/D) since thin/flat fuselages with high fineness ratios are needed to reduce drag as the vehicle moves through the atmosphere. The high-speed and long-range capabilities of waverider configurations have been demonstrated in various tests and research studies. For example, the X-51A waverider set a new record for air-breathing hypersonic flight, reaching a top speed of Mach 5.1 and flying for over six minutes [12]. Another example is the HiFire project, a joint research effort by the United States and Australia, which has conducted multiple flight tests of waverider vehicles at hypersonic speeds [13]. In such a configuration, it is natural to integrate the gas generator with the scram combustor in a top-down (or over-under) arrangement, unlike the tandem scheme sketched in Fig. 1c. The solid propellant grain in the gas generator is usually of the end-burning type with a constant burn area at all times. In the arrangement in Fig. 2, the grain usually burns from the front to the rear (with the consequent shift in vehicle center of gravity), and gaseous fuel laden with heated solid particles is injected from the fore end of the gas generator into the scram combustor. For the present study, we shall consider such a DRSJ configuration with Hydroxyl-terminated polybutadiene (HTPB) as the fuel, and a limited amount of Ammonium perchlorate (AP) as the oxidizer, with solid-particle additives like Aluminium, Magnesium, and Boron.

A typical mission would see the DRSJ vehicle launched from the ground or from an airplane in flight and boosted to a high supersonic Mach number using a solid rocket booster. The subsequent mission profile where the DRSJ engine is operational is of interest to us. We shall consider a mission profile with three segments, as marked in Fig. 3—an accelerated climb segment (AB) starting at Mach 5 at an altitude of 25 km, followed by a cruise segment (BC) at Mach 6, 30 km altitude until the fuel is exhausted at point C. The vehicle then glides without power (segment CD) for a further distance until one of the following conditions is reached: either Altitude falls to 20 km, or the vehicle decelerates to Mach 5. The net range is taken to be the total horizontal distance between points A and D. The vehicle is sized to approximately achieve a net range of 1000 km.

This typical mission profile highlights one of the key challenges in designing a working ducted rocket scramjet engine. Typically, the accelerated climb starting at the lower altitude at point A requires a much larger value of thrust, hence a higher fuel flow rate from the gas generator to the scram combustor, than the cruise segment at a higher altitude. The fuel mass flow rate into the scram combustor depends on the propellant burn rate in the gas generator and the setting
of the gas generator nozzle area. The propellant burn rate itself depends on the chamber pressure in the gas generator and the propellant burn area, which is assumed constant in this arrangement. Hence, the gas generator chamber pressure is the one variable that can be adjusted to regulate the thrust. Too high a chamber pressure can cause excessive fuel flow causing the scram combustor to experience thermal choking; on the other hand, too low a pressure may generate too little fuel flow resulting in flame out —either way, the engine fails to produce the desired thrust and the mission is jeopardized. If the system is tailored to provide a respectable rate of climb and acceleration at point A, then inevitably the minimum thrust at the 30 km cruise altitude will exceed the drag causing the vehicle to accelerate. In this process, the fuel consumption rate is larger all through, and the endpoint C where fuel is exhausted may be reached at a lower range. On the other hand, if the vehicle is designed to provide just enough thrust at 30 km altitude to balance the drag and cruise at a steady Mach number of 6, then it will not be able to generate enough thrust to get the desired climb performance at the lower Mach 5, 25 km altitude state. Thus, some sort of optimization is necessary to balance out these contradictory pulls —a solution where a compromise rate of climb is available at point A, but which also does not cause excessive acceleration between points B and C, and maximizes the range in the process. Such an exercise is formulated and carried out in the present paper.

The key parameter that needs to be controlled and monitored in this design is the gas generator chamber pressure. The only way to control the chamber pressure is by adjusting the gas generator nozzle valve area; closing the valve initially reduces the mass flow exiting the nozzle causing mass accumulation in the chamber and increased gas density and pressure. However, the higher pressure raises the propellant burn rate and increases the nozzle exit mass flow rate to more than the original value. The system attains an equilibrium at a higher value of the mass flow rates, both generated by propellant burn and exhausted from the nozzle valve. However, the relationship between the nozzle valve area setting and the chamber pressure is nonlinear and, in fact, non-minimum phase [14]. This makes the calculation of the chamber pressure that needs to be set to obtain a certain fuel flow rate somewhat complicated.

The various parameters involved in this system are intricately connected and it is not immediately obvious as to how they may be selected to obtain the ideal solution. The choice of climb angle during phase AB and the Mach number at the end of climb at point B is likely to impact the net range achieved. They are related to the thrust required for accelerated climb, which demands a certain variation of the fuel flow rate (or equivalence ratio), subject to thermal choking limits. The fuel flow rate can be obtained by regulating the gas generator chamber pressure, which, in turn, depends on the limits of the nozzle valve area variation. This also determines the minimum value of chamber pressure during cruise, which sets the rate at which the fuel is consumed and the excess of thrust over drag at cruise. As the vehicle accelerates in cruise, the lift tends to increase, and in order to maintain level flight it is necessary to compensate by lowering the angle of attack. Besides, as the fuel is consumed, the vehicle weight decreases, and so does the required lift, which calls for another downward adjustment in the angle of attack. Consequently, the lift coefficient at cruise is no longer at the optimal value to maximize range. On the other hand, if the climb angle is too shallow, the vehicle
accelerates slowly and spends too much time at lower altitudes where air density is higher, hence higher fuel flow rates are needed to maintain the same fuel equivalence ratio, causing faster depletion of the fuel. Also, the lower acceleration implies lower velocities, which in turn reduces the range covered.

It is clear from the above discussion that optimization of the various parameters in the conceptual design of a ducted rocket scramjet is not straightforward. It is necessary to model the flow and thermal processes through various components of the DRSJ to capture the effects of the parameters and the connections between them. At the same time, the model must be of reasonably low order so that it can be turned around rapidly for repeated runs during the optimization process. Computational fluid dynamics-based models can give an accurate description of the fluid-thermal processes in the engine, but usually yield models of very high order that are not useful for mission analysis or controller synthesis. On the contrary, low-order empirical models, seemingly appropriate for initial design and analysis, tend to inaccurately represent crucial physical aspects. An appropriate model is one that uses as low a mathematical complexity as possible while accurately reflecting as much of the physics as possible. Such models have been developed by us in the past for liquid-fuel ramjet engines [15, 16], and a liquid-fueled dual-combustor ramjet [17]. In these studies, models for each component were individually developed, followed by their integration to form a comprehensive engine model. This method offers the advantage of achieving high-fidelity modeling for each component while enabling the creation of a unified engine model that maintains a low complexity. In this work, a similar exercise is carried out for the first time to construct a model of the ducted rocket scramjet engine. The engine model is coupled into a flight mechanics model and used to determine the optimum value of the system parameters that maximizes the range while satisfying the various constraints imposed. The results of the optimization exercise are evaluated by conducting a time-simulated run of the mission profile for confirmation and to look for possible discrepancies.

The remainder of this paper is organized as follows: Section II provides a detailed explanation of the analytical models for each scramjet engine component. Section III focuses on the dynamics of the gas generator, describing how it maintains the required pressure via the throat valve as the solid fuel burns. In Section IV, the flight mechanics model is introduced and integrated with the scramjet engine components to create a single model for the scramjet vehicle. The problem formulation and solution methodology are also presented. Section V discusses testing the integrated model with optimized parameters using Simulink, providing real-time insights into the mission profile. Lastly, Section V analyzes the results, including external factors such as altitude, climb angle, angle of attack, and vehicle weight, etc. along with internal properties like gas generator performance and engine flow.

## II. Scramjet engine component models

The dimensions of the scramjet vehicle have been adopted from Ref. [11], and are as shown in Fig. 4, where the various DRSJ station numbers are also marked. Station ‘i’ marks the free stream conditions, while Station 1 corresponds to the vehicle nose, also the first intake ramp. The intake is a mixed-compression type comprising of two ramps of
<figure><img src="image_7.png" /><figcaption>Fig. 4 Dimensions of the scramjet vehicle with station numbers indicated.</figcaption></figure>

angles 5 deg and 10 deg, respectively. The external compression process occurs between Stations 1 and 2. Station 2 represents the start of the isolator, which extends up to the inlet of the combustor as a constant-area duct. Station 3 designates the isolator exit, which is also the entry point to the scram combustor. Station 4 signifies the exit of the combustor and the start of the aft-body nozzle, while Station ‘e’ denotes the exit plane of the nozzle. The combustor itself consists of a constant area section followed by an aft diverging section with the upper wall angled at 11 deg. The entire aft end of the vehicle functions as a single (one-sided) expansion ramp nozzle (SERN), providing a substantial exit-to-throat area ratio. The lower boundary of the nozzle flow is given by the free surface formed by the interaction of the scramjet exhaust with the hypersonic airflow around the vehicle. Additional geometric, inertial, and aerodynamic data for the configuration in Fig. 4 are provided in Appendix A.

The SERN nozzle in Fig. 4 has a ramp of 11 deg (same as the combustor aft section), which is estimated to be insufficient to expand the flow adequately to generate the necessary thrust. Instead, we shall replace the straight-line ramp of 11 deg with a curved surface which better expands the flow. This is described in a later section of this paper.

In the following, the steady-state flow model for each component of the ducted rocket scramjet engine —intake, isolator, fuel injection, scram combustor, and nozzle —is detailed. These are used for the range optimization exercise.

## A. Intake

The external compression over the two intake ramps is modeled using inviscid shock-expansion theory. Knowing the ramp angle $\theta$ and the upstream Mach number $M_u$, the shock angle $\beta$ at a ramp may be calculated using the relation:

$$ \tan \theta = 2 \cot \beta \frac{M_u^2 \sin^2 \beta - 1}{M_u^2 (\gamma_{air} + \cos 2\beta) + 2} \quad (1) $$
Subsequently, properties downstream of a shock—Mach number ($M_d$), temperature ($T_d$), density ($\rho_d$)—and the static ($P_d/P_u$) and total ($P_{0d}/P_{0u}$) pressure ratios across a shock may be obtained using the relations:

$$M_d = \sqrt{\frac{1 + \frac{\gamma-1}{2} M_u^2 \sin^2 \beta}{(\gamma M_u^2 \sin^2 \beta - \frac{\gamma-1}{2}) \sin^2(\beta - \theta)}} \quad (2)$$

$$T_d = T_u \frac{[2\gamma_{air} M_u^2 \sin^2 \beta - (\gamma_{air} - 1)][(\gamma_{air} - 1)M_u^2 \sin^2 \beta + 2]}{(\gamma_{air} + 1)^2 M_u^2 \sin^2 \beta} \quad (3)$$

$$\rho_d = \rho_u \frac{(\gamma_{air} + 1) M_u^2 \sin^2 \beta}{(\gamma_{air} - 1) M_u^2 \sin^2 \beta + 2} \quad (4)$$

$$\frac{P_d}{P_u} = \frac{2\gamma_{air} M_u^2 \sin^2 \beta - (\gamma_{air} - 1)}{\gamma_{air} + 1} \quad (5)$$

$$\frac{P_{0d}}{P_{0u}} = \left(\frac{\rho_d}{\rho_u}\right)^{\left(\frac{\gamma_{air}}{\gamma_{air}-1}\right)} \left(\frac{\gamma_{air} + 1}{2\gamma_{air} M_u^2 \sin^2 \beta - (\gamma_{air} - 1)}\right)^{\left(\frac{1}{\gamma_{air}-1}\right)} \quad (6)$$

Here, the subscript 'd' implies 'downstream of the shock' and the subscript 'u' implies 'upstream of the shock'. Most of the total pressure loss occurs across the shocks which is accounted for in Eq. 6; the viscous losses are negligible.

## B. Isolator

Further compression of the intake flow takes place in the isolator through a series of oblique shock reflections. The higher back-pressure at the exit of the isolator forces the flow to separate in the isolator duct, constricting the core flow area. Hence, the viscous losses are important here, in addition to the shock losses. The isolator flow is modeled by the classical quasi-one-dimensional gas dynamics equations (see, for example, Ref. [18]) with the addition of Ortwerth's diffuser model [19]. These are a set of ordinary differential equations for duct flows with specified stream-wise distributions of area, wall shear stress, and heat release. They allow for the calculation of the change in Mach number, core flow cross-sectional area ($A_c$), and the static pressure along the length of the duct. The isolator equations, in the manner of Ref. [20], are as follows:

$$\frac{dM^2}{M^2} = -\left(1 + \frac{\gamma_{air} - 1}{2} M^2\right) \left[ \left(\frac{2}{\gamma_{air} M^2} \frac{A}{A_c}\right) \frac{dP}{P} + \left(\frac{A}{A_c}\right) \frac{4C_f dx}{D_H} + \frac{dT_t}{T_t} \right] \quad (7)$$

$$\frac{d(A_c/A)}{A_c/A} = \left[ \frac{1 - M^2 (1 - \gamma_{air} (1 - A_c/A))}{\gamma_{air} M^2 A_c/A} \right] \frac{dP}{P} + \left( \frac{1 + (\gamma_{air} - 1) M^2}{2 A_c/A} \right) 4C_f \frac{dx}{D} + \left( 1 + \frac{\gamma_{air} - 1}{2} M^2 \right) \frac{dT_t}{T_t} \quad (8)$$

$$\frac{dP}{dx} = \frac{k}{D_H} C_{f0} \left(\frac{\rho V^2}{2}\right) \quad (9)$$

where $D_H$ is the hydraulic diameter, defined as the ratio of four times the area of the duct to its perimeter; $k = 89$ is an empirical constant determined through correlation against experimental data [21]; $C_{f0} = 0.002$ is the minimum skin
friction coefficient at the initial separation point; and $\gamma_{air}$ is the ratio of specific heats of air. Equations 7-9 are solved by space marching along the length of the isolator.

### C. Fuel injection

The mass flow rate of gaseous fuel (volume occupied by the solid particles is assumed to be negligible) injected from the gas generator into the scram combustor is calculated at the gas generator nozzle throat, assuming the throat Mach number to be sonic:

$$ \dot{m}_f^* = \sqrt{\frac{\gamma_g P_c}{R_g T_c} A_t \left(\frac{2}{\gamma_g + 1}\right)^{\frac{\gamma_g+1}{2(\gamma_g-1)}}} \quad (10) $$

where $\gamma_g$ is the ratio of specific heats for the gaseous product in the gas generator and $R_g$ is the specific gas constant. The nozzle throat area setting required to obtain a given fuel mass flow rate can be derived by inverting Eq. 10. The choked nozzle velocity may be calculated using,

$$ V_f = \sqrt{\gamma_g R_g T_f} = \sqrt{\frac{2\gamma_g R_g T_c}{\gamma_g + 1}} \quad (11) $$

### D. Scram combustor

For the scram combustor model, we are primarily interested in accurately estimating two quantities: i) the heat release due to combustion of the gaseous fuel (with solid particles) being injected from the gas generator, and ii) the flow properties at the exit station of the combustor subject to the condition that thermal choking does not occur in the scram combustor. A zero-dimensional model that accounts for friction, heat release due to combustion, heat loss at the walls, and the expansion in the combustor duct area is deemed appropriate for this purpose. The equations for a zero-dimensional analysis of compressible flow in a duct involve the conservation of mass, momentum, and energy, along with the equation of state. Following the approach in Ref. [22], these equations are expressed as follows:

$$ \rho_4 U_4 A_4 = (1 + f) \rho_3 U_3 A_3 \quad (12) $$

$$ P_4 A_4 + \rho_4 U_4^2 A_4 = P_3 A_3 + \rho_3 U_3^2 A_3 + P_3 (A_4 - A_3) - \left(\frac{1}{2}\right) \rho_3 U_3^2 C_f A_w \quad (13) $$

$$ \frac{\gamma_4 R_4}{\gamma_4 - 1} (T_4 - T_3) + \left(\frac{1}{2}\right) (U_4^2 - U_3^2) = \frac{\dot{Q}_c + \dot{Q}_w}{(1 + f) \rho_3 U_3 A_3} \quad (14) $$

$$ P_4 = \rho_4 R_4 T_4; \quad P_3 = \rho_3 R_3 T_3 \quad (15) $$

where '4' refers to the combustor exit station, and '3' denotes the combustor entry. Defining the following non-dimensional variables:
$$r = \frac{\rho_4}{\rho_3}; \quad u = \frac{U_4}{U_3}; \quad AR = \frac{A_4}{A_3}; \quad \bar{A}_w = \frac{A_w}{A_3}; \quad \bar{R} = \frac{R_4}{R_3}; \quad \bar{\gamma} = \frac{\gamma_4}{\gamma_3}$$

$$\bar{Q}_c = \frac{\dot{Q}_c}{\rho_3 U_3^3 A_3}; \quad \bar{Q}_w = \frac{\dot{Q}_w}{\rho_3 U_3^3 A_3}; \quad M_3^2 = \frac{U_3^2}{\gamma_3 R_3 T_3}; \quad M_4^2 = \frac{U_4^2}{\gamma_4 R_4 T_4}$$

Equations 12-14 can be represented as follows, where Eq. 15 is incorporated:

$$r.u.AR = 1 + f \quad (16)$$

$$r.u^2.AR \cdot \left( \frac{1}{\gamma_4 M_4^2} + 1 \right) = \left( AR \cdot \frac{1}{\gamma_3 M_3^2} + 1 \right) - \left( \frac{1}{2} C_f \bar{A}_w \right) \quad (17)$$

$$\frac{1}{\gamma_2 - 1} \left( u^2 \cdot \frac{1}{M_4^2} - \bar{R} \cdot \bar{\gamma} \cdot \frac{1}{M_3^2} \right) + \left( \frac{1}{2} \right) (u^2 - 1) = \frac{\bar{Q}_c + \bar{Q}_w}{1 + f} \quad (18)$$

Equations 16-18 are a set of nonlinear algebraic equations in three variables —r, u, M₄—which may be solved to obtain the density, velocity, and Mach number at the combustor exit (station 4). Pressure is determined using Eq. 15, while the temperature is calculated based on the Mach number.

The flow properties of the fuel-air mixture in the scram combustor are calculated by an averaging process as follows, where subscripts 'a' and 'f' refer to the airflow from the intake and the gaseous fuel injected from the gas generator, respectively:

$$\begin{gathered} \dot{m} = \dot{m}_a + \dot{m}_f, \qquad V = \frac{\dot{m}_a V_a + \dot{m}_f V_f}{\dot{m}}, \qquad \gamma = \frac{\dot{m}_a \gamma_a + \dot{m}_f \gamma_f}{\dot{m}}, \\ R = \frac{\dot{m}_a R_a + \dot{m}_f R_f}{\dot{m}}, \qquad T_0 = \frac{\dot{m}_a C_{pa} T_{0a} + \dot{m}_f C_{pf} T_{0f}}{\dot{m} C_p} \end{gathered} \quad (19)$$

The remaining quantities may then be obtained as below:

$$\rho = \frac{\dot{m}}{AV}, \quad C_p = \frac{\gamma R}{\gamma - 1}, \quad \frac{T_0}{T} = 1 + \frac{\gamma - 1}{2} M^2, \quad p = \rho RT \quad (20)$$

where M = V/a and a = $\sqrt{\gamma_g RT}$.

The heat release due to combustion of the gaseous fuel, including the solid particles, $\dot{Q}_c$, is determined using Cantera [23], as detailed below. Cantera is used as it integrates well with the rest of the models in MATLAB®. The wall heat loss $\dot{Q}_w$ is estimated to be 7.5% of heat release $\dot{Q}_c$ as per Ref. [24].

## 1. Heat release calculation

The fuel used in the gas generator is a mixture of Ammonium Perchlorate (40%), HTPB (25%), Aluminum (13%), Magnesium (12%) and Boron (10%) by mass. The pyrolysis products formed in the gas generator have been listed in
<figure><img src="image_8.png" /><figcaption>Fig. 5 Effect of temperature on Q at different equivalence ratio. The initial conditions T = 700 K, 900 K and 1100 K and P = 0.4 MPa</figcaption></figure>

Appendix B. Of these, the top 14 species by mass fraction, comprising 99.6% of the mass of the pyrolysis products are used in the heat release calculation. Of that, 60.5% consists of gaseous species and the remaining 39.5% by mass is made up of condensed species. The heat release due to the gaseous fuel is computed using Cantera, which does not handle condensed species. The additional contribution from the condensed species is evaluated using CEA (Chemical Equilibrium with Applications) code<fnref n="2" /> <fnref n="3" />, which is then used to correct the gaseous-only estimate obtained from Cantera.

A sample result from Cantera for the heat release (Q) over a range of equivalence ratios from 0.2 to 1.4, with initial pressure of 0.4 MPa, and for three values of initial temperature (700, 900, 1100 K), from the gaseous species alone, is plotted in Fig. 5. The heat release may be seen to range roughly between 9-10 and 18 MJ/kg of fuel, depending on the value of the equivalence ratio. At a typical value of equivalence ratio of 0.8, the heat release is around 12-14 MJ/kg (without considering the condensed species) depending on the scram combustor temperature at the point of fuel injection. Similar data are generated for varying values of pressure (0.6 and 0.8 MPa) for fixed values of temperature; however, very little variation is seen in the heat release with initial pressure conditions.

To correct for the effect of the condensed species, the same calculation is repeated using CEA (Chemical Equilibrium with Applications) code considering all the top 14 species by mass fraction (gaseous as well as condensed species). The results are compared with those obtained previously from Cantera for the gaseous species alone after scaling the gaseous fuel mass appropriately so that the heat release in MJ per kg of fuel is analogous in the two cases. The additional heat release ($\Delta Q_{rel}$) due to inclusion of condensed species, for a range of initial temperature and pressure conditions, is plotted in Fig. 6 for different values of the fuel equivalence ratio. Roughly, an additional heat release of 1 to 2.5 MJ/kg of fuel is obtained due to the combustion of the condensed species in the scram combustor, depending on the value of the equivalence ratio. For a particular value of the equivalence ratio equal to 0.8, the additional heat release is close to

<footnote n="2">²https://www1.grc.nasa.gov/research-and-engineering/ceaweb/</footnote>

<footnote n="3">³http://akrmys.com/public/cea/cea_install_cmd.html.en</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 6 Additional heat release ($\Delta Q_{rel} = Q_{rel,Fuel1F} - Q_{rel,Fuel1G}$) due to addition of condensed species</figcaption></figure>

2.5 MJ/kg, fairly independent of the initial temperature and pressure conditions in the scram combustor.

## 2. Check for thermal choking

When the maximum energy deposition in the burner is reached, the pressure rise in the combustion chamber due to heat release also reaches a limit. This occurs when the supersonic flow becomes critical, a condition known as thermal choking. In thermal choking, further heat addition increases the temperature, which raises the local speed of sound. When the local speed of sound exceeds the local velocity, the flow shifts toward subsonic, causing deceleration and a reduction in mass flow rate. As the mass flow rate decreases, pressure builds up, potentially leading to inlet *unstart*. Therefore, to prevent inlet unstart and maintain a supersonic Mach number at the combustor exit, it is essential to avoid thermal choking.

The scramjet flow model represented by Eqs. 16-18 may be used to detect the condition for thermal choking in the combustor duct, as explained in Ref. [akn,sreesankar arXiv]. To check for thermal choking at the combustor exit station 4, for instance, we set the Mach number there, $M_4$, as equal to 1.0. The heat release variable, $\bar{Q}_c$, is now specified as the third unknown in Eqs. 16-18, along with $r$ and $u$, and may be solved for under the condition that $M_4 = 1$ at the station with area $A_4$. Thus, a critical value of the heat release, $\bar{Q}_c^*$ is obtained corresponding to the thermally choked state. As long as the heat release value calculated from the combustion process does not exceed the critical value, the combustor is safe from thermal choking. Heat release values in excess of the critical value are untenable, and those domains in parameter space may be excluded from further analysis.

## E. Nozzle

The flow exiting the combustor expands through the Single Expansion Ramp Nozzle (SERN), as schematically depicted in Fig. 7. At the upper wall, expansion waves are generated at the corner of the single ramp. Another expansion
<figure><img src="image_10.png" /><figcaption>Fig. 7 Schematic of the expansion process through Single Expansion Ramp Nozzle (SERN)</figcaption></figure>

takes place at the tip of the nozzle cowl where the wall boundary layers along the cowl surfaces are shed in the form of a shear layer. The shear layer forms the virtual lower wall of the nozzle, separating the hot nozzle gas flow from the free stream flow around the external scramjet vehicle. The ramp expansion waves reflect off the shear layer, and vice versa, as indicated by the dashed lines in Fig. 7.

As mentioned earlier, the nozzle ramp angle provided in the original design in Fig. 4 is likely too small to provide adequate expansion. Hence, to ensure that sufficient thrust is available, the nozzle is redesigned with a curved surface replacing the straight-line ramp in Fig. 4. A method of characteristics (MOC) code [25] is used to design the curved nozzle wall and to calculate the nozzle flow properties. A minimum length nozzle is first designed with the flow expanding from Mach 1 at the nozzle throat to Mach 3.2 at the theoretical exit (not Station 'e' in Fig. 4). The shape of the ideal nozzle is shown in Fig. 8a as the upper curve at which the various characteristic lines terminate. To get a practical nozzle length, the curved shape in Fig. 8a is terminated at $y = 0.6$ m, which corresponds to the height of the rearmost point on the aft-body of the scramjet vehicle in Fig. 4. The truncated nozzle shape is marked by the bold curve in Fig. 8a. Thus, the curved nozzle fits exactly between the starting and end points of the original straight-line ramp of 11 deg as shown in 8b.

## 1. Thrust calculation

After obtaining the nozzle exit properties, the thrust produced by the nozzle flow is computed as follows:

$$T = \dot{m}_a ((1 + f)U_e - U_a) + (P_e - P_a)A_e \quad (21)$$

Equation 21 requires the exit area $A_e$ to be calculated, which is the area bounded by the upper solid wall and the lower shear layer at station 'e'. This requires the shape of the shear layer to be predicted for each operating condition of the
scramjet engine, which calls for a certain effort in cost and time that is not affordable since the model needs a quick turn-around time when used iteratively during optimization. Instead, the lower nozzle wall is assumed to be a mirror image of the upper solid wall; that is, the nozzle geometry is assumed to be axisymmetric. The error due to replacing the shear layer by an axisymmetric lower surface has been estimated for one typical flow condition —the following conditions at the combustor exit station '4' were used:

$$M_4 = 1.6, \quad P_4 = 31000 \text{ Pa}, \quad T_t = 1392.5 \text{ K}, \quad \gamma_4 = 1.2693$$

The values of $\gamma_4$ and $R_4$ are derived from Cantera [23] during the heat release calculation. Results from a two-dimensional numerical simulation which captured the shear layer were compared with those from a method of characteristics computation which assumed the lower surface to be an axisymmetric solid wall. The error in the average Mach number at the exit station 'e' was found to be 3.73%, which is quite acceptable. Hence, all further computations of the nozzle flow are carried out with the method of characteristics code, assuming the nozzle to be axisymmetric and the lower surface to be a solid wall.

In practice, the actual thrust is often lower than that calculated using Eq. 21. Usually, a thrust correction factor ($\zeta_F$) is employed, which is defined as the product of a velocity correction factor ($\zeta_v$) and a discharge correction factor ($\zeta_d$):

$$\zeta_F = \zeta_v \times \zeta_d \qquad (22)$$

These are estimated using the method described in Ref. [26]. The velocity correction factor ($\zeta_v$) is defined as the ratio

<figure><img src="image_11.png" /><figcaption>Fig. 8 Curved nozzle ramp designed using MOC to replace the 11 deg SERN ramp in Fig. 4: (a) Characteristic lines depicting the theoretical minimum length nozzle, and (b) Truncated nozzle contour obtained after terminating at y = 0.6 m (Station 'e')</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 9 Control volume analysis of gas generator</figcaption></figure>

of the kinetic energy of the exhaust leaving the nozzle to the kinetic energy of a theoretical ideal exhaust. The standard value for $\zeta_v$ is 0.92. The discharge correction factor ($\zeta_d$) is the ratio of the mass flow rate exiting from the nozzle to the ideal mass flow rate, expressed as

$$ \zeta_d = \frac{\dot{m}_{act}}{\dot{m}_{id}} \qquad (23) $$

A value of $\zeta_d = 1.01$ is considered for the present study. Therefore, the corrected thrust equation is:

$$ T = (\dot{m}_a ((1 + f)U_e - U_a) + (P_e - P_a)A_e) \times \zeta_F \qquad (24) $$

Finally, to determine the axial component of the thrust force, the axial component of exit velocity is extracted from the MOC computation. This is done by averaging the axial velocity at each point of the MOC computation at the nozzle exit plane. Likewise, an average value of the pressure is also obtained, and these are used in Eq. 24. Therefore, the axial component of the thrust, that is of interest, is calculated as follows:

$$ T_{ax} = \left( \dot{m}_a \left( (1 + f) \frac{\sum_i U_{ei} \cos(\theta_i)}{i} - U_a \right) + \left( \sum_i \frac{P_{ei}}{i} - P_a \right) A_e \right) \times \zeta_F \qquad (25) $$

where $i$ is each point on the nozzle exit plane where the MOC code returns a value, $\theta_i$ is the flow direction at each of these points, and $P_{ei}$ is the corresponding pressure.

### III. Gas generator dynamics model

The key parameter in the gas generator (GG) model is the chamber pressure $P_c$ and its time variation. To derive the relevant equation, we consider a control volume (CV) as marked in Fig. 9 which spans the end-burning grain surface, the walls of the gas generator, and the exit nozzle valve up to the throat. The rate of change of gas mass in the CV can be written in terms of the fluxes at the boundaries of the CV:

$$ \dot{m}_c = (\dot{m}_{gi} + \dot{m}_{si}) - (\dot{m}_{ge} + \dot{m}_{se}) \qquad (26) $$
where the subscripts 'g' and 's' stand for the gas and solid phase, respectively, and subscripts 'i' and 'e' denote the injection of mass at the burning surface and the mass exiting at the nozzle throat, respectively. By differentiating the relation for the perfect gas law, an alternative expression for the rate of change of mass in the GG chamber can be written as:

$$ \dot{m}_c = \frac{V_c}{R_c T_c} \frac{dP_c}{dt} + \rho_c \frac{dV_c}{dt} \qquad (27) $$

where $V_c$, $P_c$, and $T_c$ are volume, static pressure, and temperature in the GG chamber, respectively. It is reasonable to assume that the solid particles do not accumulate in the CV; that is, the mass flow rate of the solid particles injected into the CV at the burning surface is at all times equal to their mass flow exiting the CV. With this proviso, equating Eqs. 26 and 27 yields:

$$ \frac{dP_c}{dt} = \frac{R_c T_c}{V_c} \left( (\dot{m}_{gi} - \dot{m}_{ge}) - \rho_c \frac{dV_c}{dt} \right) \qquad (28) $$

The gaseous mass injected due to propellant burn is given by:

$$ \dot{m}_{gi} = \rho_b A_b \dot{r}_{burn} \qquad (29) $$

where the burn rate is modeled according to Vielle's law as $\dot{r}_{burn} = a P_c^n$, with $a = 5.1e^{-6}$ m/s, $n = 0.453$, and $P_c$ to be used in Pa. The choking mass flow rate exiting the nozzle at the throat is as follows:

$$ \dot{m}_{ge} = P_{c0} A_t \sqrt{\frac{\gamma_c}{R_c T_{c0}}} \left(\frac{2}{\gamma_c + 1}\right)^{\frac{\gamma_c + 1}{2(\gamma_c - 1)}} \qquad (30) $$

where $P_{c0}$ and $T_{c0}$ are the stagnation pressure and temperature, respectively, and $\gamma_c$ is the ratio of specific heats of the gaseous products in the gas generator. It is fair to assume that the Mach number of the flow in the gas generator is low; that is, one can assume $M_c \approx 0$. Hence, the stagnation pressure and temperature in Eq. 31 may be replaced with the corresponding static values, $P_c$ and $T_c$, respectively, as below:

$$ \dot{m}_{ge} = P_c A_t \sqrt{\frac{\gamma_c}{R_c T_c}} \left(\frac{2}{\gamma_c + 1}\right)^{\frac{\gamma_c + 1}{2(\gamma_c - 1)}} \qquad (31) $$

As the propellant in the gas generator gets burned, its quantity decreases, which may be represented by the reduction in its length:

$$ \frac{dl_{prop}}{dt} = -\dot{r}_{burn} \qquad (32) $$

Further, the GG chamber volume increases with time as the fuel is consumed, which is modeled by:

$$ \frac{dV_c}{dt} = A_b \dot{r}_{burn} \qquad (33) $$
Now, using Eqs. 31, 29, and 33 to fill in for the terms on the right-hand side of Eq. 28, the final form of the dynamical equation for the GG chamber pressure may be obtained:

$$ \frac{dP_c}{dt} = \frac{R_c T_c}{V_c} \left( -P_c A_t C_1 + C_2 P_c^n - C_3 P_c^{n+1} \right) \quad (34) $$

where

$$ C_1 = \sqrt{\frac{\gamma_c}{R_c T_c}} \left( \frac{2}{\gamma_c + 1} \right)^{\frac{\gamma_c + 1}{2(\gamma_c - 1)}}, \quad C_2 = \rho_b A_b a, \quad C_3 = \frac{A_b a}{R_c T_c} $$

are constants for the given gas generator configuration. Notably, Eq. 34 is a nonlinear differential equation for the time variation of $P_c$.

### 1. Inverse relation for valve throat area

In general, when the gas generator is operational, the chamber pressure will vary as the chamber volume is constantly increasing according to Eq. 33. If the chamber pressure is to be held constant, the nozzle valve throat area must be continuously adjusted. So an updated value of the throat area $A_t$ needs to be computed at every instant of the operation. This requires a *dynamic* inversion of the nonlinear Eq. 34. To this end, following the method in Ref. [27], the time rate of change of $P_c$ is modeled as follows:

$$ \left. \frac{dP_c}{dt} \right|_{\text{mod}} = \omega_p (P_{c\text{SET}} - P_c) \quad (35) $$

where the symbol 'mod' indicates that it is a modeled value, the subscript 'SET' refers to the required value of $P_c$, and $\omega_p$ (1/sec) is a constant, often called the bandwidth. For this study, $\omega_p = 1000$ Hz has been selected. By using Eq. 35 on the left-hand side of Eq. 34, and inverting Eq. 34 to solve for $A_t$, the required value of the throat area may be determined which will maintain the chamber pressure at the set value $P_{c\text{SET}}$:

$$ A_t = \left( \left. \frac{-V_c}{R_c T_c} \frac{dP_c}{dt} \right|_{\text{mod}} + C_2 P_c^n - C_3 P_c^{n+1} \right) \frac{1}{C_1 P_c} \quad (36) $$

## IV. Range optimization

The objective of this conceptual design exercise is twofold: 1) to confirm that the ducted rocket scramjet vehicle can indeed be operated under accelerated climb as well as level cruise conditions, using the same gas generator, and 2) to determine the optimal flight and engine parameters for such a mission (refer Fig. 3) that will yield maximum range for the given vehicle configuration and quantity of onboard fuel. For this purpose, we shall first list the flight mechanics equations to be used in this study, followed by the formulation and solution of the range optimization problem.
## A. Flight mechanics model

The flight mechanics equations are written in terms of the wind-axis variables for total velocity and flight path angle $(V, \gamma)$, following the standard practice [28]:

$$ \frac{dV}{dt} = \frac{g}{W} \left( T \cos \alpha - \frac{1}{2} \rho V^2 S C_D \right) - \frac{g}{V} \frac{dh}{dt} \qquad (37) $$

$$ \frac{d\gamma}{dt} = \frac{g}{WV} \left( \frac{1}{2} \rho V^2 S C_L + T \sin \alpha - W \cos \gamma \right) \qquad (38) $$

where the thrust force is assumed to act along the body X-axis. Both $C_L$ and $C_D$ are modeled as cubic polynomial functions of angle of attack $\alpha$, following the data in Appendix A, for Mach numbers 5, 6, and 7, and linearly interpolated between these Mach number values. The equations for range and climb rate follow:

$$ \frac{dx}{dt} = V, \qquad \frac{dh}{dt} = V \sin \gamma \qquad (39) $$

Additionally, the vehicle weight is being reduced at a variable rate depending on the GG chamber pressure, as follows:

$$ \frac{dW}{dt} = -\rho_b A_b g \times a P_c^n, \qquad a = 5.1e^{-6}, \qquad n = 0.453 \qquad (40) $$

### 1. Inverse relation for angle of attack

For the accelerated climb segment (Phase AB in Fig. 3), the aim is to gain velocity while climbing at a constant flight path angle. For the cruise segment BC, one would like to maintain a fixed altitude by holding the flight path angle uniformly to zero. However, during these flight phases, as the fuel is burned, the vehicle weight falls, and hence the lift-weight balance is affected. In order to maintain the set value of flight path angle $\gamma$, the angle of attack needs to be continuously adjusted so that Eq. 38 is balanced at all times. This requires a dynamic inversion of Eq. 38. Once again, following the method in Ref. [27], the time rate of change of $\gamma$ is modeled as follows:

$$ \left. \frac{d\gamma}{dt} \right|_{\text{mod}} = \omega_g (\gamma_{\text{SET}} - \gamma) \qquad (41) $$

where the symbol 'mod' indicates that it is a modeled value, the subscript 'SET' refers to the required value of $\gamma$, and $\omega_g$ (1/sec) is a constant bandwidth that is set to 10 Hz for this study. By using Eq. 41 on the left-hand side of Eq. 38, one obtains the function $f$ that needs to be set to zero so that $\gamma$ is driven to its set value at every instant:

$$ f(\alpha) = 0.5 \rho V^2 S C_L + T \sin(\alpha) - W \left( \left. \frac{V}{g} \frac{d\gamma}{dt} \right|_{\text{mod}} + \cos(\gamma) \right) \qquad (42) $$
The required value of angle of attack α is determined as a root of this function f(α).

## 2. Integrating flight mechanics and engine models

Consider any arbitrary point on the flight trajectory along the mission profile in Fig. 3 with the present value of Mach number, Altitude and the other flight and engine parameters known. Let the desired flight path angle at that point be γ<sub>SET</sub>, zero during the cruise segment BC and otherwise during the climb AB —this is the first Design Parameter. Solving Eq. 42 gives the required value of angle of attack at that instant. We assume that any value of angle of attack in the range [-2, 10] deg is attainable by appropriate use of the control surfaces. The vehicle drag can also then be calculated. Passing over to the engine model now, with the values of Mach number, Altitude, and Angle of attack, the intake model in Sec.II.A gives the Mach number, air mass flow rate, and other flow properties at the start of the isolator. Note that these calculations depend on the prevailing value of angle of attack; as the angle of attack varies in flight, the output of the intake model also varies accordingly. Next, the isolator model in Sec. II.B is used to obtain the flow properties at the entry to the scram combustor. The second Design Parameter is the gas generator chamber pressure P<sub>c</sub>, which may lie in the range [1.4, 7.5] MPa between the lower limit for fuel burnout and the higher limit for chamber over-pressure. The nozzle valve throat area required to maintain this value of P<sub>c</sub> may be derived from Eq. 36. With this information, the fuel injection parameters, namely, fuel mass flow rate and injection velocity, may be obtained from the model in Sec. II.C. Thereafter, the flow properties of the air-fuel mixture in the scram combustor follow from Eq. 19. Also, the fuel Equivalence Ratio (ER) may be calculated, as indicated in Appendix B, and using the tabulated data there, one can also read off the heat release rate from combustion Q̇<sub>c</sub>. Ideally, the Equivalence Ratio should lie in the range [0.4, 1.2]. The scram combustor model in Sec. II.D then allows one to check for possible occurrence of thermal choking. If thermal choking is detected, then that value of P<sub>c</sub> (and any higher value for that combination of flight and engine parameters) is inadmissible. Otherwise, the scram combustor model calculates the flow conditions at the exit of the combustor, which are then fed to the nozzle model in Sec. II.E. The nozzle model supplies the calculated value of thrust using Eq. 24.

Excess of thrust over drag allows the vehicle to accelerate or climb, or both. This is made obvious by re-framing Eq. 37 as follows:

$$ \frac{V}{W} \left( T \cos \alpha - \frac{1}{2} \rho V^2 S C_D \right) = \frac{V}{g} \frac{dV}{dt} + \frac{dh}{dt} \qquad (43) $$

where the term on the left-hand side is called the Specific Excess Power, and the two terms on the right-hand side represent the vehicle acceleration and climb rate, respectively.

Finally, during the glide descent phase CD, all the fuel has been consumed and engine thrust is zero. Maximum range during glide is obtained for a value of angle of attack corresponding to maximum lift-over-drag, as outlined in Ref. [29]. Hence, the lift-to-drag ratio (C<sub>L</sub>/C<sub>D</sub>) is to be maximized at every instant, and the corresponding value of
angle of attack used in Eqs. 37 and 38 to estimate the velocity and glide angle.

## B. Problem formulation

The immediate objective is to maximize the net range over the two mission legs, climb (AB) and cruise (BC), in Fig. 3. The two design parameters at hand are the (constant) climb angle $\gamma$ and the gas generator chamber pressure $P_c$. In practice, these two parameters are adjusted in flight by use of the angle of attack $\alpha$ via Eq. 42, and through the GG nozzle valve throat area $A_t$ as given by Eq. 36, respectively. However, for the optimization exercise, the precise values of the angle of attack $\alpha$ and the GG nozzle valve throat area $A_t$ are not relevant, as long as they remain between their respective upper and lower bounds.

The strategy to maximize the net range can be understood by considering the fuel consumed during the climb phase (AB) as a surrogate objective. If the climb angle parameter is too small, then the vehicle is in a shallow climb spending more time at lower altitudes where the air density is relatively higher, hence the fuel flow requirement at the same Equivalence Ratio condition is larger. In such a case, all the fuel may be consumed in the climb without the vehicle reaching point B. Then the cruise leg (BC) is ruled out and the net range is likely to be smaller. For an intermediate value of the climb angle, if the GG chamber pressure $P_c$ is too high, then one of two eventualities may be expected. Either the fuel injected in the scramjet combustor is too much, causing thermal choking and making the combustor inoperable, or all the fuel is consumed in accelerating the vehicle in climb leaving none for the cruise, as before. On the other hand, for a set value of GG chamber pressure $P_c$, if the climb angle parameter is selected to be too large, the specific excess power (SEP) in Eq. 43 is mostly used in gaining altitude and the vehicle acceleration is relatively smaller. Thus, the Mach number at the start of cruise (Point B) may be too small so that, along with the decreased air density at the cruise altitude, the air mass flow supplied to the scramjet combustor could be insufficient to maintain cruise flight. Therefore, one can expect a limited range of values of the two parameters, climb angle and GG chamber pressure, for which valid solutions are available.

## C. Solution methodology

A solution methodology that considers the variation of the parameters along the climb leg (AB) turns out to be too complex to be implemented during a conceptual design study. During the climb, as the air density falls, for a given equivalence ratio, the fuel flow requirement decreases, which means that the GG chamber pressure must be constantly reduced. Additionally, as the angle of attack is regulated to maintain the constant set value of the climb angle, the vehicle drag changes affecting the thrust-drag balance, and hence the vehicle acceleration. Rather than account for these effects at each evaluation of the objective function, a conservative approximation is adopted by using the value at the start of climb (Point A) as applicable for the entire climb segment AB. As a result, the estimate of fuel consumed during climb is much larger than actual as the higher value of fuel flow rate at Point A is utilized at all points in the
<table><caption>**Table 1 Optimum parameters**</caption><thead><tr><th></th><th>Climb (AB)</th><th>Cruise (BC)</th></tr></thead><tbody><tr><td>Climb angle</td><td>1.25 deg</td><td>0</td></tr><tr><td>P<sub>c</sub></td><td>3.02 MPa</td><td>1.4 MPa</td></tr><tr><td>Angle of attack at start of climb/ cruise</td><td>1.5 deg</td><td>2.6 deg</td></tr><tr><td>Optimum Mach number at B/C</td><td>5.55</td><td>5.95</td></tr><tr><td>Maximum range</td><td>243.25 km</td><td>241.43 km</td></tr></tbody></table>

climb segment. Consequently, the range evaluated by this approximation is expected to be smaller than the actual range. Nevertheless, the optimal values of the design parameters obtained by this approximation are expected to be close (if not identical) to the actual values. The true value of the optimal range for this selection of values of the design parameters is evaluated numerically in the next section, along with the actual time variation in all the system variables during the entire mission profile.

At this stage, the range is calculated by scanning the design parameter space for various combinations of values of climb angle $\gamma$ and GG chamber pressure $P_c$. The climb angle is taken to vary between 0.5 and 5 deg, while $P_c$ at the start of climb (at Point A) is varied between 1.4 and 7.5 MPa. For the cruise leg, Eqs. 37 and 38 are solved with $\gamma = 0$ and $P_c$ set at its lower limit, $P_c = 1.4$ MPa. This inevitably results in the vehicle accelerating from $M_B$ to a value around Mach 6 at point C.

The optimization calculations are carried out using a MATLAB® script along the lines outlined above. The results are shown in Fig. 10 where three parameters of interest have been plotted—the net range covering the climb and cruise segments, the Mach number at the start of cruise (Point B), and the fuel consumed by the end of climb (Point B). The constraints that were anticipated earlier are clearly visible in these plots—they are shown by the hashed regions in Fig. 10. Below a certain value of $P_c$ at Point A, there is simply not enough thrust to climb at the selected value of climb angle—this is marked by the dashed line in Fig. 10—hence, no solutions are possible. For values of $P_c$ beyond this margin, one can notice from Fig. 10c that for low values of climb angle, the entire fuel is consumed during climb and the vehicle never attains the cruise altitude of 30 km—these regions in parameter space are also discarded. For higher values of climb angles, a lesser percentage of the fuel is used up during climb, however, this improvement is cut short due to two factors. At higher values of $P_c$, the combustor encounters thermal choking making it inoperable. On the other hand, for higher values of climb angle, as seen in Fig. 10b, the Mach number at start of cruise (Point B) progressively reduces, which also decreases the air mass flow entering the intake at cruise. Below a certain value of Mach number, seen to be around 5.55, there is not enough air mass flow available to the scram combustor for operation. Consequently, as seen in Fig. 10a, the best range is obtained at a point where both the constraints—maximum $P_c$ and the minimum allowable Mach number at B—operate simultaneously. The values of the various parameters at this optimal solution are tabulated in Table 1. The optimal climb angle turns out to be 1.25 deg with $P_c$ at start of climb
<figure><img src="image_13.png" /><figcaption>Fig. 10 Results of optimization study shown as contour plots in γ – P<sub>c</sub> space, (a) net range in climb plus cruise, (b) Mach number at point B, and (c) percentage of fuel consumed by point B.</figcaption></figure>

(Point A) equal to 3.02 MPa. The corresponding value of Mach number at start of cruise is 5.55. At cruise, P<sub>c</sub> is set to the minimum allowable value of 1.4 MPa, which still leads to accelerated flight. The Mach number at end of cruise (Point C) where all the fuel has been consumed turns out to be 5.95, slightly short of the desired cruise Mach number of 6. A net range of approximately 485 km is achieved, which, however, is an underestimate as explained earlier. In practice, lesser fuel will be consumed during climb, leaving more fuel for a longer cruise phase, hence, one may expect the cruise range as well as the Mach number at end of cruise (Point C) to increase.
<figure><img src="image_14.png" /><figcaption>Fig. 11 Schematic of integrated DRSJ simulation model with the individual engine and flight mechanics blocks linked.</figcaption></figure>

## V. Numerical simulation of optimal mission profile

The prediction from the previous section is verified by conducting a numerical simulation of the integrated flight mechanics and engine model for the chosen values of the optimal parameters in Table 1. Climb angle is fixed as 1.25 deg and $P_c$ at start of climb (Point A) is set to 3.02 MPa, and all other parameters are calculated at every instant of the mission profile from the simulation. The various sub-models of the ducted rocket scramjet, as detailed in Secs. II and III, are suitably linked to create a single, unified model for the DRSJ engine in Simulink®. This model is then integrated with the flight mechanics model, as described in Sec. IV.A, and a final integrated Simulink® model is created, as shown in Fig. 11. The flight mechanics model is placed within the block labeled 'Variation in freestream properties.' Atmospheric properties are obtained by the COESA atmospheric model, and $g$ is determined using the WGS84 exact gravity model available in Simulink®. A variable-step stiff solver ode23s, based on a modified Rosenbrock formula of order two, is used.

The intake mass flow rate is corrected for flow spillage at the cowl under off-design conditions by a multiplicative flow factor. A check for the possible occurrence of thermal choking at the scram combustor exit is placed —if the exit Mach number $M_4$ falls below 1.1, then the set value of the gas generator chamber pressure $P_{cSET}$ is decreased in steps
<figure><img src="image_15.png" /><figcaption>Fig. 12 Variation in flight mechanics parameters with time</figcaption></figure>

of 0.5 MPa until the thermal choking condition is relieved. This causes the GG valve throat area to be appropriately adjusted, reducing the fuel flow rate into the combustor, and hence the Equivalence Ratio, which should result in a fall in the heat release rate in the scram combustor. The GG chamber pressure is also lower limited at 1.4 MPa to avoid flame burnout in the gas generator. The powered phase terminates at point C in Fig. 3 when the propellant length diminishes to 3% of its initial size, which is taken to be the unburned sliver fraction. Thereafter, the GG properties (except the sliver propellant length) are artificially set to zero, and the glide phase is initiated. A stopping criterion is placed which halts the simulation when either the altitude in glide falls to 20 km or the vehicle slows down to Mach 5.

## A. External properties

Figure 12 shows the time variation of the altitude, climb angle, vehicle velocity, and weight over the entire mission (climb, cruise, and glide). The DRSJ commences its climb from an altitude of 25 km, achieving 30 km in roughly 2 minutes. The vehicle operates at 30 km until the total consumption of fuel, which occurs approximately 6 minutes into the flight. Subsequently, the vehicle glides to an altitude of around 29 km, coinciding with achieving a flight Mach of 5. The glide phase lasts for approximately 2 minutes. For the entire climb phase, a constant climb angle of 1.25deg (optimum value) is maintained, which is reset to 0 at cruise to hold a constant altitude. During the glide, the descent angle is observed to increase with time as the AOA is constantly adjusted to ensure that the lift-to-drag ($C_L/C_D$) ratio is maximized at each instant. The vehicle continues to accelerate through the climb and cruise phases. It can be observed that the glide phase starts at around 320 seconds, by which time the weight is the empty weight of the vehicle without
fuel.

Figure 13 represents the time variation of the angle of attack (AOA), L/D, intake air mass flow rate, and thrust. During climb and cruise, the AOA is calculated by the inverse relationship in Eq. 42 so as to maintain the respective climb angle. As altitude increases, the air mass flow rate reduces due to decreasing density. This decrease in density surpasses the impact of the rising velocity causing the air mass flow entering the intake to fall steadily during the climb. In level flight, where density remains constant, the air mass flow rate escalates in tandem with the increasing velocity. The thrust is observed to decrease as the vehicle climbs up, primarily due to the effect of lower air density. During level flight, with constant $P_c$, the thrust produced exceeds drag, causing an increase in velocity (Fig. 12).

## B. Internal properties

The variation in some of the gas generator (GG) properties with time is illustrated in Fig. 14. The GG chamber pressure ($P_c$) is initially set at 3.02 MPa (optimum $P_c$ for climb). During the climb, there are several instances where a decrease in the set value of $P_c$ is noted, indicating that the scram combustor likely registered a Mach number near 1, indicating thermal choking. Otherwise, the set $P_c$ during the climb is maintained by scheduling the valve throat area as per Eq. 36—this is also noted in Fig. 14. For level flight, $P_c$ is set to the minimum value of 1.4 MPa, to minimize the acceleration during cruise. The fuel flow rate adjusts itself in tandem with $P_c$, and the Equivalence Ratio (ER) varies between 0.4 and 0.6 depending on the mutual variation of the fuel and air mass flow rates.

<figure><img src="image_16.png" /><figcaption>Fig. 13 Variation of vehicle velocity, intake air mass flow rate, weight, and L/D with time</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 14 Plots depicting variations of some of the GG properties with time</figcaption></figure>

Fig. 15 represents plots of Mach number vs time at the freestream ($M_i$), intake exit ($M_2$), isolator exit ($M_3$), scram combustor exit ($M_4$), and nozzle exit ($M_e$) stations. The vehicle is seen to accelerate in climb, and then further accelerate during cruise up to Mach 6.2 before it runs out of fuel. The flow is slowed down in the intake, and the Mach number is further reduced to 2.3-2.5 by the end of the isolator (entry to scram combustor). At the combustor exit, $M_4$ approaches sonic several times —those are the instances when the set $P_c$ is reduced to prevent thermal choking. The Mach number at the start of cruise (point B) from Fig. 15 is quite close to 5.55 —the value estimated during the optimization exercise in the previous section, whereas $M_C = 6.2$, slightly larger than the 5.95 estimate in Table 1.

The integrated range from the simulation is plotted in Fig. 16—the DRSJ covers approximately 230 km during the climb, 314 km in level flight, and adds another 232 km to its range during the glide phase. The total range achieved is approximately 776 km. Interestingly, the climb+cruise range from Fig. 16 is 544 km, which is not very significantly different from the 485 km estimate from the optimization study. This suggests that the approximations considered for the optimization solution were not unreasonable.
<figure><img src="image_18.png" /><figcaption>Fig. 15 Plot of Mach number vs time</figcaption></figure>

VI. Conclusion

A fundamental issue in the conceptual design of a solid-fuel ducted rocket scramjet (DRSJ) vehicle is the contrary requirements of high thrust at start of climb and relatively lower thrust during higher-altitude cruise. In addition, there are constraints on the available intake air mass flow rate, and avoidance of thermal choking in the scram combustor. The challenge of optimizing the mission range of a solid-fuel DRSJ under these requirements and constraints has been undertaken in this work. Appropriate models for the various components of the scramjet engine, along with equations for the flight trajectory, and the pressure dynamics in the gas generator, have been integrated to simulate the DRSJ operation. Optimum values of the design parameters —climb angle and GG chamber pressure at start of climb —have
<figure><img src="image_19.png" /><figcaption>Fig. 16 Plot of range vs time</figcaption></figure>

been computed, and the optimum is found to lie at the boundary defined by the thermal choking limit and the intake air mass flow rate limit at start of cruise. The optimal solution is verified by a time simulation of the mission profile with the angle of attack and the GG nozzle throat area continuously varied to maintain the optimum values of climb angle and GG chamber pressure throughout the flight. The integrated aerothermal-flight mechanics model introduced in this work is key to the design and optimization of such hypersonic vehicles.

# Appendices

## A. Geometric, inertial, and aerodynamic data for the scramjet vehicle

The geometric and inertial data for the DRSJ vehicle, adapted from Ref. [11], are listed in Table 2. These are classified under four headings: Airframe, Engine, Wings, and Weights. Additionally, aerodynamic data in the form of lift coefficient, $C_L$, and drag coefficient, $C_D$, as a function of angle of attack for free stream Mach number values of 5, 6, and 7, are reproduced from Ref. [11] in Fig. 17. It is to be noted that these aerodynamic estimates are obtained in Ref. [11] using inviscid flow analysis. Therefore, a multiplicative viscous correction factor of 1.092 based on our internal validation studies is applied uniformly for all Mach numbers to the $C_D$ data in Fig. 17.

## B. Fuel species mass fractions

The pyrolysis products of the fuel in the gas generator of the DRSJ are listed in Table 3. Species mass and mole fractions are provided and they are arranged in decreasing order of their mole fractions.
<table><caption>Table 2 Geometric data for scramjet vehicle</caption><thead><tr><th>System</th><th>Variable</th><th>Value</th><th>Meaning</th></tr></thead><tbody><tr><td rowspan="6">Airframe</td><td>L<sub>body</sub></td><td>4.27 m</td><td>Body length</td></tr><tr><td>W<sub>body</sub></td><td>0.5842 m</td><td>Body width</td></tr><tr><td>H<sub>upsurface</sub></td><td>0.2 m</td><td>Body height</td></tr><tr><td>W<sub>engine</sub></td><td>0.4 m</td><td>Engine width</td></tr><tr><td>y<sub>engine</sub></td><td>-0.4 m</td><td>Longitudinal coordinate of engine bottom</td></tr><tr><td>S<sub>ref</sub></td><td>0.2986 m<sup>2</sup></td><td>Reference area</td></tr><tr><td rowspan="6">Engine</td><td>δ1<sub>forebody</sub></td><td>5°</td><td>First wedge inclination of forebody</td></tr><tr><td>δ2<sub>forebody</sub></td><td>10°</td><td>Second wedge inclination of forebody</td></tr><tr><td>x1<sub>forebody</sub></td><td>0.2</td><td>Axial relative coordinate of intersection of wedge 1 and 2</td></tr><tr><td>x2<sub>forebody</sub></td><td>0.4</td><td>Axial relative coordinate of intersection of wedge 2 and scramjet engine</td></tr><tr><td>x<sub>nozzle</sub></td><td>0.75</td><td>Axial relative coordinate of intersection of scramjet engine and afterbody</td></tr><tr><td>x1<sub>engine</sub></td><td>0.4</td><td>Axial relative coordinate of front point of engine bottom</td></tr><tr><td rowspan="5">Wings</td><td>x2<sub>engine</sub></td><td>0.75</td><td>Axial relative coordinate of rear point of engine bottom</td></tr><tr><td>L<sub>wingroot</sub></td><td>0.94 m</td><td>Chord length of wing root</td></tr><tr><td>L<sub>wingtip</sub></td><td>0.5589 m</td><td>Chord length of wing tip</td></tr><tr><td>L<sub>wingspan</sub></td><td>0.33 m</td><td>Wing span</td></tr><tr><td>θ<sub>sweepback</sub></td><td>45°</td><td>Sweep angle of leading edge</td></tr><tr><td rowspan="2">Weights</td><td>c<sub>wingthickness</sub></td><td>0.02 m</td><td>Relative thickness</td></tr><tr><td>Vehicle</td><td>671.33 kg</td><td>Vehicle weight with fuel</td></tr><tr><td></td><td>Fuel</td><td>226.02 kg</td><td>Fuel weight</td></tr></tbody></table>

References

[1] Xiang, Z., Zhixun, X., Likun, M., Chaolong, L., Chuanbo, F., Natan, B., and Alon, G., "Research progress on solid-fueled Scramjet," Chinese Journal of Aeronautics, Vol. 35, No. 1, 2022, pp. 398-415.

[2] Vaught, C., Witt, M., Netzer, D., and Gany, A., "Investigation of solid-fuel, dual-mode combustion ramjets," Journal of Propulsion and Power, Vol. 8, No. 5, 1992, pp. 1004-1011.

[3] Lv, Z., Xia, Z.-x., Liu, B., and Liu, Y.-C., "Experimental and numerical investigation of a solid-fuel rocket scramjet combustor," Journal of Propulsion and Power, Vol. 32, No. 2, 2016, pp. 273-278.
<table><caption>Table 3 Pyrolysis products formed in gas generator with molar and mass fractions</caption><thead><tr><th>Species</th><th>Mass fraction</th><th>Mole fraction</th><th>Species</th><th>Mass fraction</th><th>Mole fraction</th></tr></thead><tbody><tr><td>H2</td><td>0.040799</td><td>0.403430</td><td>ALH</td><td>0.000253</td><td>0.000180</td></tr><tr><td>CO</td><td>0.258610</td><td>0.184050</td><td>ALOH</td><td>0.000287</td><td>0.000130</td></tr><tr><td>C(gr)</td><td>0.091275</td><td>0.151490</td><td>ALCL2</td><td>0.000491</td><td>0.000100</td></tr><tr><td>Mg</td><td>0.086553</td><td>0.070991</td><td>CH4</td><td>0.000072</td><td>0.000090</td></tr><tr><td>BN(cr)</td><td>0.082376</td><td>0.066171</td><td>BO</td><td>0.000067</td><td>0.000050</td></tr><tr><td>ALCL</td><td>0.166740</td><td>0.053241</td><td>BOCL</td><td>0.000094</td><td>0.000030</td></tr><tr><td>B4C(cr)</td><td>0.075415</td><td>0.027210</td><td>ALCL3</td><td>0.000134</td><td>0.000020</td></tr><tr><td>MgAL2O4(cr)</td><td>0.144160</td><td>0.020200</td><td>ALHCL2</td><td>0.000099</td><td>0.000020</td></tr><tr><td>HBO</td><td>0.012042</td><td>0.008630</td><td>HBO2</td><td>0.000044</td><td>0.000020</td></tr><tr><td>MgCL2</td><td>0.021730</td><td>0.004550</td><td>HNC</td><td>0.000027</td><td>0.000020</td></tr><tr><td>MgCL</td><td>0.007943</td><td>0.002650</td><td>MgH</td><td>0.000025</td><td>0.000020</td></tr><tr><td>HCL</td><td>0.004609</td><td>0.002520</td><td>ALHCL</td><td>0.000032</td><td>0.000010</td></tr><tr><td>H</td><td>0.000059</td><td>0.001170</td><td>BH2CL</td><td>0.000024</td><td>0.000010</td></tr><tr><td>AL2O</td><td>0.002527</td><td>0.000720</td><td>BH3</td><td>0.000007</td><td>0.000010</td></tr><tr><td>N2</td><td>0.000956</td><td>0.000680</td><td>CL</td><td>0.000018</td><td>0.000010</td></tr><tr><td>AL</td><td>0.000677</td><td>0.000500</td><td>H2O</td><td>0.000009</td><td>0.000010</td></tr><tr><td>C2H2_acetylene</td><td>0.000536</td><td>0.000410</td><td></td><td></td><td></td></tr><tr><td>B2O2</td><td>0.000888</td><td>0.000330</td><td></td><td></td><td></td></tr><tr><td>HCN</td><td>0.000434</td><td>0.000320</td><td></td><td></td><td></td></tr></tbody></table>

[4] Lv, Z., Xia, Z.-x., Liu, B., and Huang, L.-y., "Preliminary experimental study on solid-fuel rocket scramjet combustor," *Journal of Zhejiang University-SCIENCE A*, Vol. 2, No. 18, 2017, pp. 106-112.

[5] Huang, L., Hu, G., Hu, B., Tian, L., and Yang, Y., "Experiment on combustion of solid rocket scramjet," *Journal of Solid Rocket Technology*, Vol. 43, No. 200, 2020, pp. 5-9.

[6] Liu, J., Wang, N.-f., Wang, J., and Li, Z.-y., "Optimizing combustion performance in a solid rocket scramjet engine," *Aerospace Science and Technology*, Vol. 99, 2020, p. 105560.

<figure><img src="image_20.png" /><figcaption>Fig. 17 Variation of C<sub>L</sub> and C<sub>D</sub> with angle of attack, for Ma = 5, 6, and 7 [11]</figcaption></figure>
[7] Li, C., Xia, Z., Ma, L., Zhao, X., and Chen, B., "Numerical study on the solid fuel rocket scramjet combustor with cavity," *Energies*, Vol. 12, No. 7, 2019, p. 1235.

[8] Li, C., Zhao, X., Xia, Z., Ma, L., and Chen, B., "Influence of the vortex generator on the performance of solid rocket scramjet combustor," *Acta Astronautica*, Vol. 164, 2019, pp. 174-183.

[9] Bauer, C., Davenne, F., Hopfe, N., and Kurth, G., "Modelling of a throttleable ducted rocket propulsion system," *47th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit*, 2011, p. 5610.

[10] Rathi, N., and Ramakrishna, P., "Attaining hypersonic flight with aluminum-based fuel-rich propellant," *Journal of Propulsion and Power*, Vol. 33, No. 5, 2017, pp. 1207-1217.

[11] Li, H., Ping, L., and Dajun, X., "Control-oriented modeling for air-breathing hypersonic vehicle using parameterized configuration approach," *Chinese Journal of Aeronautics*, Vol. 24, No. 1, 2011, pp. 81-89.

[12] Musielak, D., *Scramjet Propulsion: A Practical Introduction*, John Wiley & Sons, 2022.

[13] Kimmel, R., Adamczak, D., Berger, K., and Choudhari, M., "HIFiRE-5 flight vehicle design," *40th Fluid Dynamics Conference and Exhibit*, 2010, p. 4985.

[14] Alan, A., Yildiz, Y., and Poyraz, U., "Gas generator pressure control in throttleable ducted rockets: A classical and adaptive control approach," *51st AIAA/SAE/ASEE Joint Propulsion Conference*, 2015, p. 4236.

[15] Gupta, N. K., Gupta, B. K., Ananthkrishnan, N., Shevare, G. R., Park, I. S., and Yoon, H. G., "Integrated modeling and simulation of an air-breathing combustion system dynamics," *AIAA modeling and simulation technologies conference and exhibit*, 2007, p. 6374.

[16] Chandra, K. P. B., Gupta, N. K., Ananthkrishnan, N., Park, I. S., and Yoon, H. G., "Modeling, simulation, and controller design for an air-breathing combustion system," *Journal of Propulsion and Power*, Vol. 26, No. 3, 2010, pp. 562-574.

[17] Wadwankar, N., Kandasamy, G., Ananthkrishnan, N., Renganathan, V., Park, I.-S., and Hwang, K.-Y., "Dual combustor ramjet engine dynamics modeling and simulation for design analysis," *Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering*, Vol. 233, No. 4, 2019, pp. 1307-1322.

[18] Shapiro, A. H., *The dynamics and thermodynamics of compressible fluid flow*, New York: Ronald Press, 1953.

[19] Ortwerth, P., "Scramjet vehicle integration," *Scramjet Propulsion*, 2001, pp. 1105-1289.

[20] Smart, M. K., "Flow modeling of pseudoshocks in backpressured ducts," *AIAA Journal*, Vol. 53, No. 12, 2015, pp. 3577-3588.

[21] Smart, M. K., "Scramjet isolators," *AVT-185 AVT/VKI Lecture Series, von Kármán Inst. RTO-EN-AVT-185-AC/323 (AVT-185) TP/375, Rhode St. Genèse, Belgium*, 2010.

[22] Wu, C., Zhong, F., and Fan, J., "An analytical theory of heated duct flows in supersonic combustors," *Theoretical and Applied Mechanics Letters*, Vol. 4, No. 3, 2014, p. 032001.
[23] Goodwin, D. G., Moffat, H. K., and Speth, R. L., “Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes,”, 2018.

[24] O’Brien, T. F., Starkey, R. P., and Lewis, M. J., “Quasi-one-dimensional high-speed engine model with finite-rate chemistry,” *Journal of Propulsion and Power*, Vol. 17, No. 6, 2001, pp. 1366–1374.

[25] Zucrow, M. J., and Hoffman, J. D., *Gas dynamics. volume 2-multidimensional flow*, Wiley, 1977.

[26] Sutton, G. P., and Biblarz, O., *Rocket Propulsion Elements*, John Wiley & Sons, 2016.

[27] Looye, G., “Rapid prototyping using inversion-based control and object-oriented modelling,” *Nonlinear analysis and synthesis techniques for aircraft control*, 2007, pp. 147–173.

[28] Sinha, N. K., and Ananthkrishnan, N., *Advanced flight dynamics with elements of flight control*, CRC Press, 2017.

[29] Vinh, N. X., *Flight mechanics of high-performance aircraft*, Vol. 4, Cambridge university press, 1993.