A Review on Supersonic Nozzle Design and Analysis with Traditional Methods

1. Introduction

1.1. History of Rocket Propulsion

Rockets and their fundamentals have been around humankind history for a while, with their bigger developments being associated with warfare. Starting as bamboo tubes filled with gunpowder used in religious celebrities, primitive rockets are invented in the 1^st century China. “Arrows of flying fire”, resembling solid propellant rockets, are the first warfare application of rocket propulsion, used in 1232 by the Chinese army to repel the Mongol invaders.

It is during the 20th century that rocketry is heavily developed into almost nowadays technology, being encourage by the needs of World War Two and the Cold War.

Robert H. Goddard, an american engineer, designs and builds the first liquid propellant rocket. 34 rockets where launch between 1926 and 1941, getting to altitudes of up to 2.6 km and speeds as high as 885 km/h. Later, in 1945, the single and multi-stage rocket engines for both rocket and jet-assisted takeoff are developed by Goddard.

During the same time, german engineer Wernher von Braun develops the V-2, a liquid-propellant rocket missile used by the Axis Forces in World War Two. In the after math of the war, von Braun is recruited by NASA becoming a prominent figure in the developing of the Saturn V launcher.

During the height of the Space Race, driven by the Cold War, the development of rocket technology reached its peak, resulting in several notable milestones[1]:

• 1957 – The launch of the first Earth-orbiting artificial satellite, Sputnik I, by the USSR. The satellite Sputnik II carried a dog named Laika into orbit, later that year.
• 1959 – The lunar surface is firstly touched by a human-made object, the USSR’s Luna 2
• 1961 – Venus is first flew by the USSR’s Venera 1, being the first to reach another planet.
• 1962 – Yuri Gagarin is the first human in space onboard the USSR’s Vostok spacecraft.
• 1969 – The American Apollo 11, powered by Saturn V rockets, successfully landed the first humans on the moon, Neil Armstrong and Buzz Aldrin.
• 2013 – Voyager 1, launched in 1977 by the USA, becomes the first human-made object to leave the solar system.

Nowadays, several nations like Europeans countries (under ESA), India, Japan, China, Brazil, among many others, have their own space programs, for which rockets, and their development, are crucial. Even more recently, the space industry stopped being a centralized state research to become also a private business, as recent Space X, Blue Origin and Virgin Galactic launches prove.

2. De Laval Convergent-Divergent Nozzle

Thrust is generated in aerial vehicles by mass ejection, as per Newton’s third law, mathematically by the simplified Equation 1. Besides the first parcel of the Equation corresponding to mass ejection, a second force is generated by the pressure difference between the flow at the exit and the ambient pressure, translated in the second parcel [1].

$$F=\dot{m}{u}_e+(P_e-P_{amb})A_e$$

(1)

Nozzles where invented as means to change the properties of a flow such as velocity and pressure. The de Laval nozzle, containing a convergent followed by a divergent section does so in such a way that a supersonic flow can be obtain at the exit. The flow behavior off such nozzle is observed in Figure 1.

A de Laval nozzle is composed of three sections: the convergent (subsonic), the throat (sonic) and the divergent (supersonic). In each regions, the flow behaves differently, with differentiated analysis methods to determinate its contribution for the production of thrust.

The convergent section accepts the hot and higly pressured exhaust gases from the combustion chamber, at almost static conditions, providing the first mean of acceleration until a sonic value. So, the convergent affects the mass flow and, at some degree, the efficiency of the combustion chamber.

At the throat, the flow reaches a sonic velocity, being choked, meaning it can not be accelerated further without area increases. Having a fixed area and Mach Number, it’s the throat, and flow’s properties dependent upon the conditions provided by the combustion chamber, that limits the mass flow through the nozzle, posing great impact on the thrust produced according to the first parcel of Equation 1.

The same parcel of thrust’s Equation 1, is also greatly impacted by the exit velocity. Since the flow is choked at the throat, is in the divergent section that accelerates it further to supersonic values. The comprehensible flow keeps on gaining kinetic energy as a trade off by losing internal energy (as seen in Figure 1 by the decrease in temperature), with a pressure drop. The exit velocity is dependent on the ratio between exit and throat areas, outside conditions in relation to the chamber and the divergent section design that affects the flow’s orientation and can induce losses.

Three major concepts are necessary to evaluate a nozzle’s performance. The first concept, represented by Equation 2, is the specific impulse, which measures the thrust produced by the weight flow rate of propellant burned. It is important to understand the rocket’s nozzle efficiency.

$$I_{sp}={\displaystyle \frac{F}{\dot{m}{g}_0}}$$

(2)

The second concept, represented by Equation 3, is the total impulse, which measures the thrust produced over the entire duration of the engine’s operation.

$$I_t={\int }_0^{t}Fdt$$

(3)

Finally, the thrust coefficient, exhibited in Equation 4, represents the dimensionless value of thrust obtained by dividing it by the chamber pressure and the throat area. This is a fair indicator of nozzle efficiency since, in theory, the thrust value depends solely on the pressure ratio between the chamber and the ambient, and the area ratio between the exit plane and the throat [3].

$$C_F={\displaystyle \frac{F}{p_1{A}_t}}$$

(4)

2.1. Ideal Nozzle

The ideal nozzle is a concept that provides the best theoretical performance for the given conditions, to which all nozzle designs are compared to. A parallel flow and exit pressure matching external pressure are some consequences from this simplification. Looking into more detail, the ideal nozzle involves the following considerations [3]:

• The exhaust gases are in chemical equilibrium, with a homogeneous composition and in a gaseous state.
• The exhaust gases obey the perfect gases law.
• The flow is adiabatic, meaning there are no appreciable heat transfers to the wall.
• The flow is isentropic, without discontinuities in properties and/or shock waves.
• The boundary layer effects are neglected, with no wall friction decelerating parts of the flow.
• The flow rate is constant and steady, without gas pulsations, turbulence, and with the transient effects of starting-up and shutting down neglected since they are of short duration.
• The exit flow is parallel and uni-dimensional.
• The flow’s velocity, pressure, temperature, and density are uniform across any nozzle’s normal section.
• Before the combustion process, propellants are store at ambient temperature, except for cryogenics which are at their boiling point.

Nonetheless, this are fair simplifications due to how fast the expansion occurs. For chemical rockets, tested efficiency parameters are usually just 1%-6% below those calculated with totally ideal assumptions [3].

3. Real Nozzle

The level of precision required for rocket systems and to optimize an already very efficient system, demands more accurate algorithms that involve a better understanding of energy losses, heat transfer and other physical or chemical phenomena.

In real nozzles, not all the chemical energy contained in the propellant is available as thermal energy, and not all the thermal energy can be converted into kinetic energy. Therefore, the performance of a rocket nozzle is affected by several factors and non-idealities that can cause losses in thrust and specific impulse.

The first factor, already stated before as a concern, is the divergency losses resulting from a non-unidimensional flow profile at the exit plane due nozzle contours that don’t end in a null slope. Another factor is that small chamber areas relative to the throat, lead to pressure losses in the chamber and a reduction of the exhaust exit velocity and thrust produced.

Boundary layer effects can also be significant, particularly for smaller nozzles with low area ratios, making a portion of the flow, ranging from $2\%$ to $25\%$, subsonic. Due to the no-slip condition, the flow next to the wall has zero-speed with high thermal energy that is dissipated to the wall and the nearby moving flow. The immediate flow is laminar and subsonic, turning further from the wall into transonic/supersonic zones where turbulence can occur, before it joins the totally supersonic free stream.

Combustion in the rocket propulsion is typically unsteady, resulting in flow oscillations and lower chamber pressure during transient operations. Chemical reactions in the flow can cause losses of up to $0.5\%$, and incomplete mixing of the gas composition means that the gas constants are not uniform throughout the nozzle, contributing for some mathematical incongruities of the numerical models.

Small solid particles and liquid droplets also exist in the flow, which need to be accelerated, with little thermal energy contribution. If bigger than 0.015mm and equivalent to more than 6% of the flow’s mass, specific impulse losses can rise to $10\%$-$20\%$. Although less significant, in high area ratio nozzles, aggressive expansion can lead to the precipitation of some exhaust constituents.

It still must be mentioned that uncooled materials suffer erosion, a phenomena of upmost relevance for the throat since, if its area increases, some chamber pressure is lost, with an additional area ratio decrease. Up to $0.7\%$ of specific impulse can be lost.

And, at last, operating away from the design altitude reduces thrust and specific impulse for nozzles with fixed area ratio [3].

4. Nozzle Phenomena

4.1. Operation Critical Points

Fixing the chamber pressure, flow begins in subsonic regime throughout the entire nozzle when the outside pressure starts to drop, expanding in the divergent and xompressing in the divergent. As the pressure ratio ($P_t/P_e$) increases, the flow accelerates and the mass flow rate rises until it reaches a sonic value at the throat, causing the flow to choke. Beyond the throat, the flow keeps compressing in the divergent section, and it returns to a subsonic regime. This point is called the first critical point and is represented in Figure 2 by curve a. Any infinitesimal decrease in the outside pressure beyond the first critical point will cause the flow in the divergent section to change from subsonic to supersonic.

Now that there is a supersonic flow in the divergent section, it no longer compresses but expands, but the great area variation promotes a more aggressive expansion than the one needed to equal the ambient pressure. Thus, a shock forms inside the divergent section, starting at the throat and moving to the exit plane as the outside pressure drops and a non-isentropic compression besides the shock is needed so, after the shock, the flow becomes subsonic and uses the remaining divergent section as a compressor. This is represented by line b in Figure 2. When the outside pressure is low enough, the shock reaches the exit plane, as suggested by line c in the same figure, representing the second critical point.

After the second critical point, the flow will not be able to expand much further even if the ambient pressure keeps dropping. However, it still reaches the exit plane at higher pressures than the ambient. Therefore, the normal shock at the exit begins to turn oblique, with its intensity/angle reducing as the outside pressure decreases. This phenomenon is known as an overexpanding nozzle.

The oblique shock will eventually reach a $0^{\circ }$ angle, representing no shock at all, and the nozzle reaches its third critical point. At this point, the flow is isentropic and is used as the nozzle’s design point.

Further decreases in outside pressure will cause the flow to have a pressure higher than the ambient pressure at the exit, which is known as underexpansion. In this case, expansion fans will be generated from the nozzle’s lips to allow the flow to be further expanded beyond the exit plane.

4.2. Overexpansion and Underexpansion Phenomena

The overexpansion occurs when the ambient pressure is higher than the one at the third critical point. A weak shock is needed at the exit so that the flow pressure increases until it matches the ambient pressure.

But the oblique shock will deviate the flow, so a second shock is needed to realign the flow. This will increase the pressure, making it higher than the ambient pressure, so now an expansion fan is needed to expand the flow and match the ambient pressure once again.

The expansion fan also deviates the flow, requiring another expansion fan to straighten the flow, returning it to a state where its pressure is again below the ambient pressure.

The process repeats itself, as shown in Figure 3b), forming a structure known as the Mach Diamond. The underexpansion process, depicted in Figure 3a), is similar, but starts with a pair of expansion fans at the exit .

5. Theory of Flow in Nozzles

6. Introduction to High Velocity Compressible Flow

When a flow subjected to great gradients and accelerates to high velocities from stagnation, density cannot be assumed to be constant raising the concept of compressible flow. Examining Equation 5, if the pressure experienced by a finite volume of fluid increases, compression occurs, resulting in a reduction in volume and an increase in density. The symbol $\tau$ represents the compressibility of the gas, which has different values for isothermal and isentropic processes.

$$\tau =-{\displaystyle \frac{1}{vol}}{\displaystyle \frac{dvol}{dP}}$$

(5)

High-speed flows are often initiated and propelled by strong pressure gradients. Therefore, the impact of pressure on the gas density cannot be ignored. Equation 6 relates the three main properties for any two states in a gaseous flows assuming isentropy. [6].

$${\displaystyle \frac{P_1}{P_2}}={\left({\displaystyle \frac{{\rho }_1}{{\rho }_2}}\right)}^{\gamma }={\left({\displaystyle \frac{T_1}{T_2}}\right)}^{{\displaystyle \frac{\gamma }{\gamma -1}}}$$

(6)

6.1. Speed of Sound and Mach Number

All particles within a fluid move in random directions with a certain velocity, attributing the fluid a specific internal energy. This random kinetic propagation speed is known as the speed of sound, and it is dependent on the fluid and its internal energy, as described by Equation 7 for thermally and calorically perfect gases. When a perturbation is applied to the fluid, it will excite nearby molecules, which will eventually collide with other molecules, transmitting the disturbance at a specified mean speed.

$$a=\sqrt{\gamma RT}$$

(7)

The Mach number is defined as the ratio between the velocity of the flow and the velocity of sound, as shown in Equation 8. It represents the relationship between the flow’s kinetic energy and the random molecular kinetic energy. Based on it, different flow regimes can be established, including:

• Incompressible (M < 0.3) – Density variations are small, mostly because these flows are not associated with strong pressure gradients. Hence, this flow can be assumed to possess constant density.
• Subsonic (0.3 < M < 0.8 at freestream) – Property variations are always continuous, and the flow exhibits straight and parallel streamlines that move, converge and shape around any obstacle. In every point, the flow has a Mach number less than 1, but compressibility effects cannot be ignored.
• Transient (0.8 < M< 1.2) – As the flow approaches the speed of sound, it is not possible to guarantee that all the flow is subsonic since it may accelerate while contouring an object, creating supersonic "pockets." Shocks can appear, indicating discontinuous changes in properties.
• Supersonic (M > 1) – The entire flow moves faster than the speed of sound. Streamlines do not bend around objects, except when encountering a shock or experiencing an expansion wave. After a shock, the flow must remain supersonic.
• Hypersonic (M > 5) – At such high speeds, a shock wave causes explosive changes in properties. The temperature can increase so much that molecular dissociation effects must be considered. Nothing particularly special happens at M = 5, as the referred phenomenon increases with the Mach number, being just a convention.

$$M={\displaystyle \frac{V}{a}}$$

(8)

The Mach number allows the calculation of the Mach Angle, as shown in Equation 9, which is half of the Mach Cone’s angle. This structure is formed by pressure waves and appears around bodies moving faster than the speed of sound. To be noted that the outside of the Mach Cone is called the silence zone, where the presence of a moving object cannot be heard.

$$\mu =si{n}^{-1}{\displaystyle \frac{1}{M}}$$

(9)

6.2. Stagnation Properties

All motion is described in relation to a reference frame. In compressible flow, an isentropic deceleration can be implemented until the flow becomes static in relation to such reference frame. The new values of the fluid’s properties are called total or stagnation properties, identified with the suffix “0”. In a rocket’s nozzle, this reference state is the conditions felt at the combustion chamber, as velocity is so small compared to that at the nozzle, so it can be overlooked.

The static values of the properties are associated with the stagnation ones by the velocity at which the fluid is moving, as described by Equations 10, 11 and 12. As the Mach number of the flow increases by an isentropic process, these three static properties decrease.

$${\displaystyle \frac{T_0}{T}}=1+{\displaystyle \frac{\gamma -1}{2}}{M}^2$$

(10)

$${\displaystyle \frac{P_0}{P}}={\left(1+{\displaystyle \frac{\gamma -1}{2}}{M}^2\right)}^{{\displaystyle \frac{\gamma }{\gamma -1}}}$$

(11)

$${\displaystyle \frac{{\rho }_0}{\rho }}={\left(1+{\displaystyle \frac{\gamma -1}{2}}{M}^2\right)}^{{\displaystyle \frac{1}{\gamma -1}}}$$

(12)

6.3. Normal or Strong Shock

A shock represents a discontinuity in fluid properties that occurs in a very thin region, meaning there are significant temperature and pressure gradients, where viscosity and dissipative effects are strongly felt. As a result, the process is not isentropic, and the relations in Equation 6 are no longer valid. As no heat is added or removed from the flow, the process is adiabatic.

Equations 13, 14, 15, and 16 establish the properties of the flow across the normal shock through an adiabatic and irreversible process.

$$M_2^2={\displaystyle \frac{1+\frac{\gamma -1}{2}{M}_1^2}{\gamma M_1^2-\frac{\gamma -1}{2}}}$$

(13)

$${\displaystyle \frac{{\rho }_2}{{\rho }_1}}={\displaystyle \frac{(\gamma +1)M_1^2}{2+(\gamma -1)M_1^2}}$$

(14)

$${\displaystyle \frac{P_2}{P_1}}=1+{\displaystyle \frac{2\gamma }{\gamma +1}}(M_1^2-1)$$

(15)

$${\displaystyle \frac{T_2}{T_1}}=\left[1+{\displaystyle \frac{2\gamma }{\gamma +1}}(M_1^2-1)\right]\left[{\displaystyle \frac{2+(\gamma -1)M_1^2}{(\gamma +1)M_1^2}}\right]$$

(16)

Nonetheless, these equations admit that the flow before the shock can be subsonic (M<1). Therefore, the second law of thermodynamics must be taken into account through Equation 17. For an unitary Mach number, an isentropic $(s2-s1=0)$ infinitely weak normal shock is obtained, and for a subsonic Mach number, a negative variation of entropy is predicted, which is physically impossible.

$$s_2-s_1=c_{p}ln\left({\displaystyle \frac{T_2}{T_1}}\right)-Rln\left({\displaystyle \frac{P_2}{P_1}}\right)$$

(17)

Since shocks are adiabatic processes, the stagnation temperature remains the same. Nevertheless, the stagnation pressure drops, and it is related to the entropy gain as shown in Equation 18.

$$s_2-s_1=-Rln\left({\displaystyle \frac{P_{02}}{P_{01}}}\right)$$

(18)

6.4. Oblique or Normal Shock

More broadly, shocks are considered oblique, as illustrated by Figure 4, and are formed when the flow must be deflected at a certain angle, $\theta$, due to the presence of a concave surface. The streamlines are always parallel, changing direction discretely at the shock.

$$M_{n1}=M_{1}sin\beta$$

(19)

$$M_{t1}=M_{1}cos\beta$$

(20)

The flow is decomposed into two components, one normal to the shock, according to Equation 19, which is analyzed as a normal shock, and one parallel component that is not affected by the shock, as per Equation 20. Note that $\beta$ represents the shock angle that decomposes the velocity of the incoming flow, and $\theta$ represents the deviation experienced by the flow after the shock, coinciding with the concavity of the surface.

The Mach number after the oblique shock is obtained by Equation 21. The angles $\beta$ and $\theta$ are related to each other according to Equation 22, which is graphically represented in Figure 5.

$$M_2={\displaystyle \frac{M_{2n}}{sin(\beta -\theta )}}$$

(21)

$$tan\theta ={\displaystyle \frac{2}{tan\beta }}\left[{\displaystyle \frac{M_1^2{sin}^2\beta -1}{M_1^2(\gamma +cos2\beta )+2}}\right]$$

(22)

There is a maximum angle ${\theta }_{max}$ to which the flow can be deflected by a straight oblique shock. If the surface has $\theta >{\theta }_{max}$, the shock will be bent and flow detachment may be induced.

Similarly, if the deflection angle $\theta$ is fixed, as the Mach number decreases, the shock angle increases until it reaches a value where the fixed $\theta$ equals the maximum possible deflection angle ${\theta }_{max}$ for that Mach number. Beyond this point, further decreases in the Mach number also do not have a straight shock solution.

For any given deflection angle, there are two possible shock wave angles. The smaller shock angle corresponds to the weak shock solution, with the flow after the shock remaining supersonic $(M_2>1)$. The highest shock angle corresponds to the strong shock solution, with the flow after the shock being subsonic $(M_2<1)$. The occurrence of each solution depends on the pressure conditions downstream in relation to the pre-shock flow. A strong pressure gradient will favor the strong shock solution, while a weaker gradient will favor the weak shock solution.

6.5. Prandtl-Meyer Expansion

Looking at Figure 6, when a supersonic flow encounters a convex surface, it must deflect from itself, wich promotes an expansion, as opposed to an oblique shock. Velocity increases, and the density, pressure, and temperature decrease. Also, these changes in properties are so smooth, that they can be assumed to be continuous, corresponding to an isentropic process. This occurs because the expansion is composed of a fan of waves, each of which contributes to infinitesimally accelerate the flow.

Basically, the vertex in Figure 6 emanates infinitesimal Mach lines, responsible for expanding the flow. This Mach lines are characteristic lines.

Due to the isentropic behavior, the changes in the flow’s temperature and pressure are given by Equations 23 and 24, respectively. Density can be easily calculated using Equation 6 under an ideal gas assumption.

$${\displaystyle \frac{T_2}{T_1}}={\displaystyle \frac{1+\frac{\gamma -1}{2}{M}_1^2}{1+\frac{\gamma -1}{2}{M}_2^2}}$$

(23)

$${\displaystyle \frac{P_2}{P_1}}={\left[{\displaystyle \frac{1+\frac{\gamma -1}{2}{M}_1^2}{1+\frac{\gamma -1}{2}{M}_2^2}}\right]}^{{\displaystyle \frac{\gamma }{\gamma -1}}}$$

(24)

The Prandtl-Meyer angle is a characteristic of the flow, as a function of the Mach number by the definition of Equation 25 . As for the reflection angle, it is simply given as the difference between the Prandtl-Meyer angle before and after the expansion, as stated in Equation 26.

$$\nu \left(M\right)=\sqrt{{\displaystyle \frac{\gamma +1}{\gamma -1}}}{tan}^{-1}\sqrt{{\displaystyle \frac{\gamma -1}{\gamma +1}}(M^2-1)}-{tan}^{-1}\sqrt{M^2-1}$$

(25)

$$\Delta \theta =\nu \left(M_2\right)-\nu \left(M_1\right)$$

(26)

When numerical algorithms are implemented to solve compressible flows, sometimes it is necessary to calculate the Mach number from its associated Prandtl-Meyer angle. Obviously there is no analytical solution for Equation 25 so some inversion methods were developed [7], one of which is described in Appendix A.

7. Quasi-1D Flow

Similar to unidimensional flow, quasi-one-dimensional flow represents an improvement as it also allows for variations in area. The assumption of smooth wall angles, prevents shock formation, ensuring isentropy. Although the geometry is two-dimensional, the fluid properties are still assumed to be dependent only on the direction of the moving flow, meaning that they are constant through any given cross-sectional area.

The flow is still considered adiabatic and isentropic. The continuity equation also considers the transversal area, as suggested in Equation 27. This leads to the relation between area and Mach number stated in equation 28.

$${\rho }_1{V}_1{A}_1={\rho }_2{V}_2{A}_2$$

(27)

$${\displaystyle \frac{dA}{A}}=(M^2-1){\displaystyle \frac{du}{u}}$$

(28)

For subsonic velocities, when section area decreases, velocity increases, and vice-versa. In a rocket nozzle, this results in the flow being accelerated by decreasing the area, forming the convergent segment. For supersonic velocities, area increases lead to velocity increases, with the opposite being true. Therefore, if the area increases, as in a rocket nozzle’s divergent, the velocity also increases. As for sonic velocities, Equation 29 suggests that a minimum area is reached, which in a rocket nozzle represents the throat, mathematically reasoning the chocking phenomena.

$${\displaystyle \frac{dA}{A}}=0$$

(29)

The area-Mach number relation, given by Equation 30, shows that the Mach number at any location in the duct is a function of the ratio between the local duct area and the sonic throat area. There are two solutions, a subsonic and a supersonic one. As this is an isentropic process, the relations for this type of process can be easily applied to solve the flow properties at any point in the duct.

$${\displaystyle \frac{A}{A^{*}}}={\displaystyle \frac{1}{M}}{\left[{\displaystyle \frac{2}{\gamma +1}}\left(1+{\displaystyle \frac{\gamma -1}{2}}{M}^2\right)\right]}^{{\displaystyle \frac{\gamma +1}{2(\gamma -1)}}}$$

(30)

8. Hypersonic Flow

Contemporary space and defense industries have a keen interest in developing hypersonic vehicles. In the past, hypersonic engineering challenges have already emerged; for instance, during the Apollo 11 mission, reentry into Earth’s atmosphere occurred at a velocity of Mach 36.

At Mach 5, although there isn’t a notable event akin to the sound barrier breakthrough at Mach 1, the flow field begins to be dominated by phenomena that are insignificant at lower velocities. As such, the study of hypersonic aerodynamics becomes a distinct subject within the discipline of compressible aerodynamics.

8.1. Initial Definition

In hypersonic flow, the shock formed on a wedge will be positioned extremely close to the surface, resulting in a very thin shock layer, as illustrated in Figure 7. The higher the incident Mach number, the stronger the shock, leading to a significant increase in fluid density downstream of the shock. As a consequence, the flow downstream of the shock can more easily "squeeze" into smaller spaces, resulting in a reduction in volume while conserving mass flow. Furthermore, if high enthalpy flow effects are considered, such as chemical reactions in the gas, the shock wave will have an angle even closer to that of the wedge.

The shock wave can interact with the boundary layer, causing it to thicken and leading to modeling complexities, especially for laminar boundary layers. At high Reynolds numbers, the shock layer can be considered inviscid, and Newtonian Theory provides a simple and straightforward method for modeling the fluid dynamics.

The analysis of supersonic flow reveals that entropy increases through a shock and is directly proportional to the shock strength. Consequently, when the flow encounters a blunt nose, different streamlines experience varying degrees of shock strength, resulting in strong entropy gradients around the nose. These gradients extend away from the shock, creating an "entropy layer" around the body. This layer interacts with the boundary layer, posing additional challenges for its analytical modeling.

Viscous dissipation occurs as kinetic energy is converted into internal energy due to the velocity profile imposed by the no-slip condition within the boundary layer. This phenomenon leads to an increase in the local static temperature of the flow, as illustrated by Figure 8, which shows the static temperature profile of a hypersonic boundary layer.

A hypersonic boundary layer is typically thick and develops rapidly due to the high temperatures involved. This increase in temperature leads to higher viscosity coefficients, thickening the boundary layer. Additionally, to maintain mass flow conservation, the thickness of the boundary layer must increase since there’s a decrease in density, keeping pressure constant across a transverse section.

In engineering applications, this thickening of the boundary layer causes a significant displacement of inviscid streamlines, making bodies appear much wider. This displacement can lead to changes in the freestream flow, which further exacerbates the growth of the boundary layer in a feedback loop process known as viscous interaction. This interaction affects calculations of lift, drag, and stability due to its impact on surface pressure distribution.

Friction within the boundary layer dissipates kinetic energy, increasing the static temperature to values where molecular vibration, dissociation, and ionization may occur. These phenomena challenge the assumption of a calorically perfect gas and is further complicated if the vehicle is coated with ablative shielding materials.

Strong shocks in hypersonic flow can significantly elevate the freestream flow temperature, influencing the aerodynamic behavior of the vehicle. Heat transfer rates, including convection heating and radiation from the hot gas in the shock layer, become crucial design parameters in engineering applications. The presence of hot plasma around a reentering vehicle can also pose challenges for communications.

While high enthalpy flows introduce additional complexities, simpler analysis can assume a calorically perfect gas. Hypersonic flow exhibits various phenomena, including thin shock layers, entropy layers, viscous interaction effects, and high temperatures, all of which must be considered in engineering design and analysis.

8.2. Wave Relations

Equations presented previously for shock relations are valid for any Mach number greater than unity, so they can also be applied to hypersonic flow. In fact, for hypersonic flow, some approximations and simplifications can be made, that are important to understand the flow’s behavior.

Most simplifications advent from the fact that every term of $M_{\infty }^2$ and higher absolves smaller terms when $M_{\infty }\to \infty$.

$${\displaystyle \frac{p_2}{p_{\infty }}}\to {\displaystyle \frac{2\gamma }{\gamma +1}}{M}_{\infty }^2{sin}^2\beta$$

(31)

$${\displaystyle \frac{{\rho }_2}{{\rho }_{\infty }}}\to {\displaystyle \frac{\gamma +1}{\gamma -1}}$$

(32)

$${\displaystyle \frac{T_2}{T_{\infty }}}\to {\displaystyle \frac{2\gamma (\gamma -1)}{{\left(\gamma +1\right)}^2}}{M}_{\infty }^2{sin}^2\beta$$

(33)

$${\displaystyle \frac{u_2}{V_{\infty }}}\to 1-{\displaystyle \frac{2{sin}^2\beta }{\gamma +1}}$$

(34)

$${\displaystyle \frac{v_2}{V_{\infty }}}\to {\displaystyle \frac{sin\left(2\beta \right)}{\gamma +1}}$$

(35)

$$C_p\to \left({\displaystyle \frac{4}{\gamma +1}}\right){sin}^2\beta$$

(36)

$$\beta \to {\displaystyle \frac{\gamma +1}{2}}\theta$$

(37)

Looking at previous Equations it can be seen that the temperature and pressure ratios can become infinitely large while density and velocity components ratios reach finite limits when $M_{\infty }\to \infty$.

8.3. Local Surface Inclination Method

Linearized supersonic flow leads to the pressure coefficient given by Equation 38. Therefore, the pressure coeficient is a linear function of the plate inclination, for given freestream conditions. In theory this analysis can be extended to supersonic values.

$$C_p={\displaystyle \frac{2\theta }{\sqrt{M_{\infty }^2-1}}}$$

(38)

From Equation 36, obtained as $M_{\infty }\to \infty$, a condition $\gamma \to 1$ yields $C_p\to 2{sin}^2\beta$. Similarly, applying the same condition to Equation 32 results in ${\rho }_2/{\rho }_{\infty }\to \infty$. By considering mass conservation, the shock angle aligns with the surface inclination, i.e., $\theta =\beta$. Thus, Equation 39 is derived, which coincides with the Newtonian Theory.

$$C_p=2{sin}^2\theta$$

(39)

It can be concluded that as $M_{\infty }\to \infty$ and $\gamma \to 1$, Newtonian Theory provides a better approximation for modeling the flow. However, real hypersonic flow may not adhere strictly to these conditions, rendering Newtonian Theory inaccurate when compared to experimental results. Nevertheless, the theory yields results accurate enough for preliminary analyses where simplicity is desired.

With Newtonian Theory, the aerodynamic behavior of a flat plate is described by Equations 40, 41, and 42.

$$c_l=2{sin}^2\alpha cos\alpha$$

(40)

$$c_d=2{sin}^3\alpha$$

(41)

$${\displaystyle \frac{L}{D}}=cot\alpha$$

(42)

The lift-to-drag ratio of the plate increases as the plate inclination $\alpha$ decreases, approaching infinity when $\alpha \to 0^{\circ }$. However, when friction effects are considered, drag will always have a finite value, leading $L/D\to 0$ as $\alpha \to 0^{\circ }$. Unlike subsonic and supersonic flow, lift is not linear for $0^{\circ }<\alpha <{15}^{\circ }$, but after, peaking at $\alpha =54.7^{\circ }$.

Experimental results of the pressure coefficient for various geometries compared with Newtonian Theory are illustrated in Figure 9. It is evident that Newtonian Theory becomes more accurate with increasing Mach number, particularly beyond $M_{\infty }>5$. Additionally, there is higher accuracy for 3D objects (cone) compared to 2D objects (wedge).

There are other surface inclination methods for hypersonic flow such as the tangent wedge, tangent cone and shock-expansion.

9. High Enthalpy Flow

Apollo 11 mission reentered the atmosphere at a velocity of 11 km/s, equivalent to Mach 32.5 at an altitude of 53 km. Analyzing the strong shock that forms at the blunt nose of the vehicle, assuming a calorically perfect gas, results in a shock layer with a static temperature of 58300 K, which is incorrect. When heating the air, before reaching such temperature, the gas gets vibrationally excited and can even chemically react or ionize, absorbing energy. In fact, the shock layer is composed by a partially ionized plasma. When the previous analysis is performed but with a chemically reacting gas, and the adiabatic index becomes a function of both pressure and temperature, static temperature is "just" 11600 K.

If atmospheric air is heated at a pressure of $1atm$, oxygen and nitrogen molecules will begin to dissociate at temperatures of $2500K$ and $4000K$, respectively. Once the temperature reaches $9000K$, all the molecules have fully dissociated, and atom ionization begins. From a calorically perfect gas perspective, these phenomena "consume" heat that does not contribute to an increase in the system’s temperature.

9.1. Microscopic Description of Gases

Quantum mechanics dictates that particles can only possess quantized values of energy, with the lowest value observed at $T=0K$, denoted by ${ϵ}_0^{\prime }$ and referred to as the zero-point energy, defining the ground state. The energy of a molecule can be measured relative to this ground state, such that ${ϵ}_j={ϵ}_j^{\prime }-{ϵ}_0$.

The energy of a molecule is the sum of four energy modes, as suggested by Equation 43. These energy modes are quantized, associated with a particular atomic behavior.

$${\epsilon }^{\prime }={\epsilon }_{\mathrm{trans}}^{\prime }+{\epsilon }_{\mathrm{rot}}^{\prime }+{\epsilon }_{\mathrm{vib}}^{\prime }+{\epsilon }_{\mathrm{el}}^{\prime }$$

(43)

9.1.1. Translational Energy

As depicted in Figure 10a), this energetic mode pertains to the translational motion of the center of mass of a molecule. It encompasses the associated kinetic translational energy, which can be further decomposed into three degrees of geometric and thermal freedom (one for each orthogonal direction).

All quantized modes of translational energy are closely spaced, giving the illusion of continuous changes in energy. Equation 44 quantifies this energetic mode, where $n_i$ are quantum integer numbers, h is the Planck’s constant and $a_i$ linear dimensions that describe the system’s volume.

$${\epsilon }_{\mathrm{trans}}^{\prime }={\displaystyle \frac{h^2}{8m}}\left({\displaystyle \frac{n_1^2}{a_1^2}}+{\displaystyle \frac{n_2^2}{a_2^2}}+{\displaystyle \frac{n_3^2}{a_3^2}}\right)$$

(44)

9.1.2. Rotational Energy

Besides translation, a molecule can rotate around the orthogonal axis around its center of mass, as seen in Figure 10b). Both kinetic rotational energy and the moment of inertia contribute to this energetic mode.

Similar to translational motion, rotational motion also exhibits three degrees of geometric and thermal freedom. However, linear molecules have negligible rotation around the z-axis, resulting in only two degrees of freedom.

The energy levels become more distant with increasing energy, indicating that it takes more energy to transition from one level to the next. It is the only energetic mode that has a total value of zero at the ground state.

Rotational energy is defined by Equation 45, where J is the rotational quantum number and I the moment of inertia of the molecule.

$${\epsilon }_{\mathrm{rot}}^{\prime }={\displaystyle \frac{h^2}{8{\pi }^{2}l}}J(J+1)$$

(45)

9.1.3. Vibrational Energy

Atoms within a molecule oscillate around their equilibrium positions. For a diatomic molecule, this motion can be visualized as a spring connecting the two atoms, as depicted in Figure 10c).

In this analogy, the vibrational mode involves kinetic energy associated with the linear motion of each atom and potential energy resulting from intramolecular forces. Despite having only one degree of geometric freedom, diatomic molecules exhibit two degrees of thermal freedom. However, polyatomic molecules are more complex, featuring multiple vibrational modes.

The spacing between vibrational energy levels is even greater than that of rotational energy levels, but this spacing decreases. Consequently, it requires less energy to transition from one vibrational level to the next compared to transitioning between rotational levels.

Vibrational energy is defined by Equation 46, where n is the vibrational quantum integer number and v the fundamental vibrational frequency of the molecule.

$${\epsilon }_{\mathrm{vib}}^{\prime }=hv\left(n+{\displaystyle \frac{1}{2}}\right)$$

(46)

9.1.4. Electronic Energy

An atom, as depicted in Figure 10d), consists of electrons orbiting the nucleus with a certain rotational kinetic energy and potential energy arising from their electromagnetic interaction.

In this context, the concepts of geometric and thermal degrees of freedom become obsolete. Additionally, it’s easy to understand why the total electronic energy is above zero at the ground state. If it were null, it would imply that the electrons had collapsed into the nucleus.

The spacing between each energy level in the electronic mode is the largest among all energetic modes, and it decreases with increasing energy.

9.1.5. Macrostate

The molecules that composed a thermodynamic system will mandatorily occupy quantized energy levels. $N_j$ is the number of molecules that inhabit the level of energy ${ϵ}_j^{\prime }$, constituting that level’s population. The total energy of the system is given by Equation 47.

$$E=\sum _j{\epsilon }_j^{\prime }{N}_j$$

(47)

Due to intermolecular collisions, the molecules within a system are always changing their energy level, meaning the population distribution is constantly changing. Each different population distribution is called a macrostate.

There’s a macrostate that is the most likely to occur at any given time. Statistical mechanics delivers the tools to predict the most likely macrostate, with it corresponding to thermodynamic equilibrium.

Within an energy level, there are some microstates, that its population $N_j$ can inhabit, without changing its macrostate. Molecular collisions also promote the alternating of microstates internally to the energy level.

Since each microstate is considered to take place with equal probability, the most likely macrostate is the one that has the highest number of microstates associated. Usually, there’s a single macrostate of the system that has considerably many more microstates than the others due to the high number of molecules that constitute the typical system for engineering applications.

9.1.6. Microstate

The angular momentum of a given molecule, being a vector quantity, can have different directions but all with the same associated energy ${ϵ}_i^{\prime }$. The angular momentum is also quantized, meaning only a few directions are allowed. For example, in Figure 11, only the three depicted orientations are allowed, without intermediate states. Analogous phenomenon occurs for vibrational and electronic modes.

What was described are distinguishable scenarios among the population of a given energy level. So, each possible variation within an energy level is called a degeneracy or statistical weight, and each is represented by $g_j$.

The microstate is defined as the distribution of the population among the degenerate states of a certain energy level.

There are two theoretical cases for calculating the population of each energy level for the most likely macrostate.

9.1.7. Boson

If a molecule has an even number of elementary particles, it will follow the Bose-Einstein statistical distribution.

Each microstate can be inhabited by an unlimited number of the $N_j$ molecules from the corresponding energy level.

Equation 48 represents the number of microstates of the $N_j$ indistinguishable molecules from energy level ${ϵ}_j^{\prime }$ that has $g_j$ degenerate states.

$${\displaystyle \frac{(N_j+g_j-1)!}{(g_j-1)!N_j!}}$$

(48)

Thermodynamic probability, W, can be defined as the sum of Equation 48 for the several energy levels j that the molecules of the system can occupy. Put simply, it measures the disorder of the system by counting the number of possible microstates within a macrostate.

Each macrostate will have different values of W, and the most probable macrostate is the one with the largest value of W. To find it, W must be maximized, a task from which Equation 49 emerges. The symbol * denotes the most likely macrostate and is used to define the population $N_j^{*}$ for each energy level ${ϵ}_j^{\prime }$, thereby defining a macrostate.

$$N_j^{*}={\displaystyle \frac{g_j}{e^{\alpha }{e}^{\beta {ϵ}_j^{\prime }}-1}}$$

(49)

9.1.8. Fermion

If a molecule has an odd number of elementary particles, it will obey the Fermi-Dirac statistical distribution.

Each microstate can only be inhabited by a single molecule at a time, so $g_j>N_j$. Equation 50 gives the number of microstates for an energy level.

$${\displaystyle \frac{g_j!}{(g_j-N_j)!N_j!}}$$

(50)

Analogously to the boson, the most likely macrostate for a fermion is given by Equation 51.

$$N_j^{*}={\displaystyle \frac{g_j}{e^{\alpha }{e}^{\beta {ϵ}_j^{\prime }}+1}}$$

(51)

9.1.9. Boltzmann Distribution

At low temperatures (like $T=5K$), the differences between bosons and fermions are evident as molecules densely occupy the lower levels of energy. But at high temperatures, there are many more energy levels being inhabited and so, each population is scarce, meaning $N_j<<g_j$.

So the population distribution is given for both bosons and fermions by Equation 52. This is called the Boltzmann distribution and is applicable for temperatures considerably above $0K$, like in most engineering applications.

$$N_j^{*}=g_j{e}^{-\alpha }{e}^{-\beta {\epsilon }_j^{\prime }}$$

(52)

A partition function is defined in Equation 53 that is function of the systems temperature and volume.

$$Q\equiv \sum _j{g}_j{e}^{-{\epsilon }_j/kT}$$

(53)

$\beta$ can be obtained by combining both classical and statistical mechanics in Equation 54, where k is the Boltzmann constant.

$$\beta ={\displaystyle \frac{1}{kT}}$$

(54)

From zero-point energetic measurements and some mathematical manipulation, $\alpha$ can be replaced, yielding Equation 55.

$$N_j^{*}=N{\displaystyle \frac{g_j{e}^{-{\epsilon }_j/kT}}{Q}}$$

(55)

So, in a system with N particles, quantum mechanics states that there are well-defined quantized energy levels ${ϵ}_j$ with a certain $g_j$ amount of degenerate states each, on which molecules will be distributed. Equation 55 precisely determines the number of molecules $N_j^{*}$ on each energy level ${ϵ}_j$ when the system reaches thermodynamic equilibrium at a certain temperature and volume.

9.1.10. Partition Function as a Function of Volume and Temperature

The partition function of a molecule is the product of the partition function for each energetic mode, as suggested in Equation 56. The volume dependence arises from the multiplication of the spatial dimension of the molecule $a_i$ in Equation 57. Temperature dependence is evident in all energetic modes.

$$Q=Q_{\mathrm{trans}}{Q}_{\mathrm{rot}}{Q}_{\mathrm{vib}}{Q}_{\mathrm{el}}$$

(56)

$$Q_{\mathrm{trans}}={\left({\displaystyle \frac{2\pi mkT}{h^2}}\right)}^{3/2}{a}_1{a}_2{a}_3$$

(57)

$$Q_{\mathrm{rot}}={\displaystyle \frac{8{\pi }^{2}IkT}{h^2}}$$

(58)

$$Q_{\mathrm{vib}}={\displaystyle \frac{1}{1-e^{-h\mathrm{v}/kT}}}$$

(59)

$$Q_{\mathrm{el}}\equiv \sum _{l=0}^{\infty }{g}_l{e}^{-{\epsilon }_1/kT}=g_o+g_1{e}^{-{\epsilon }_1/kT}+g_2{e}^{-{\epsilon }_2/kT}+\cdots$$

(60)

To note, Equation 59 is only applicable to diatomic molecules, as polyatomic molecules possess multiple and complex vibrational modes. Equation 60 takes into account considerations from spectroscopic data concerning the electronic levels ${ϵ}_i$. The value of ${ϵ}_1$ is relatively high, such that for $T\le 15000K$, a second-degree expansion provides accurate results.

9.1.11. Microscopic Thermodynamic Properties

In equilibrium, the internal energy of the system is obtained with Equation 61, measured above the zero-point. At constant volume, in a system with N molecules, it can be generalized to Equation 62. Often, internal energy per unit mass is desired, so Equation 63 can be used, where R represents the universal gas constant.

$$E=\sum _j{\epsilon }_j{N}_j^{*}$$

(61)

$$E=Nk{T}^2{\left({\displaystyle \frac{\partial lnQ}{\partial T}}\right)}_V$$

(62)

$$e=R{T}^2{\left({\displaystyle \frac{\partial lnQ}{\partial T}}\right)}_V$$

(63)

Other thermodynamic properties can be derived from statistical mechanics, as shown in Equations 64, 65, and 66 for enthalpy, entropy, and pressure, respectively.

$$h=RT+R{T}^2{\left({\displaystyle \frac{\partial lnQ}{\partial T}}\right)}_V$$

(64)

$$S=Nk\left(ln{\displaystyle \frac{Q}{N}}+1\right)+NkT{\left({\displaystyle \frac{\partial lnQ}{\partial T}}\right)}_V$$

(65)

$$p=NkT{\left({\displaystyle \frac{\partial lnQ}{\partial V}}\right)}_T$$

(66)

9.2. Thermodynamic Properties Evolution

The Equipartition Theorem states that each thermal degree of freedom contributes $1/2RT$ to the internal energy of a gas. This applies to both translational and rotational internal energy, as shown in Equations 67 and 68, respectively.

$$e_{\mathrm{trans}}={\displaystyle \frac{3}{2}}RT$$

(67)

$$e_{\mathrm{rot}}=RT$$

(68)

However, the Equipartition Theorem does not apply to the vibrational mode. This is because it was formulated using classical mechanics based on macroscopic observations, without considering microscopic and quantum mechanics. Typically, $e_{\mathrm{vib}}<RT$, and its value for diatomic molecules is given in Equation 69.

$$e_{\mathrm{vib}}={\displaystyle \frac{h\nu /kT}{e^{h\nu /kT}-1}}RT$$

(69)

The sensible energy of the system is the sum of the previous three energetic modes and the electronic mode $e_{el}$, for which a specific formula is not provided. Both the sensible energy and the specific heat at constant volume, depicted in Equation 70, are functions of temperature alone in a thermally perfect gas. This is due to the disregard of intermolecular interactions, which simplifies the equally likely microstates.

$$c_v={\displaystyle \frac{3}{2}}R+R+{\displaystyle \frac{{(h\nu /kT)}^2{e}^{h\upsilon /kT}}{{(e^{h\nu /kT}-1)}^2}}R+{\displaystyle \frac{\partial e_{\mathrm{el}}}{\partial T}}$$

(70)

Now, the definition of a calorically perfect gas also becomes more evident. For air, $c_v=(5/2)R$, $c_p=c_v+R=(7/2)R$, and $\gamma =c_p/c_v=7/5$. It is clear that this simplification only considers translational and rotational energies, which is fairly valid for $T\le 600K$, as there are little to no effects from the vibrational mode.

Looking at Figure 12, one could expect $c_v\to (7/2)R$ as $T\to \infty$, but this is not true, as at higher temperatures, effects from dissociation and ionization change $c_v$. Another aspect from the graphic is that at $T=0K$, there’s only a contribution from the translational mode and $c_v=(3/2)R$, with the initial increases in temperature responsible for increasing rotational energy from a null value.

9.3. Chemical Equilibrium

Consider the chemical reaction described by Equation 71, derivated from the generic form of Equation 72. Chemical equilibrium is achieved when the system has matured to the point where its composition (quantity of each species) is constant, and each direction of the chemical equation is occurring at the same pace.

$$AB\rightleftarrows A+B$$

(71)

$$\sum _{i=1}^n{\nu }_i^{\prime }{X}_i\underset{k_b}{\stackrel{k_f}{\rightleftarrows }}\sum _{i=1}^n{\nu }_i^{\prime \prime }{X}_i$$

(72)

Noting that $N^x$ is the number of molecules of the species x in the system and $N_j^x$ the number of molecules of the species x in the energy level ${ϵ}_j^x$, the total energy of the system is given by Equation 73.

$$E=E^{\mathrm{A}}+E^{\mathrm{B}}+E^{\mathrm{AB}}=\mathrm{const}$$

(73)

where:

$$E^{\mathrm{A}}=\sum _j{N}_j^{\mathrm{A}}{\epsilon }_j^{\prime \mathrm{A}}=\sum _j{N}_j^{\mathrm{A}}\left({\epsilon }_j^{\mathrm{A}}+{\epsilon }_o^{\mathrm{A}}\right)$$

(74)

$$E^{\mathrm{B}}=\sum _j{N}_j^{\mathrm{B}}{\epsilon }_j^{\prime \mathrm{B}}=\sum _j{N}_j^{\mathrm{B}}({\epsilon }_j^{\mathrm{B}}+{\epsilon }_o^{\mathrm{B}})$$

(75)

$$E^{\mathrm{AB}}=\sum _j{N}_j^{\mathrm{AB}}{\epsilon }_j^{\prime \mathrm{AB}}=\sum _j{N}_j^{\mathrm{AB}}\left({\epsilon }_j^{\mathrm{AB}}+{\epsilon }_o^{\mathrm{AB}}\right)$$

(76)

The microscopic calculation of chemical equilibrium can be translated, once again, to finding the most likely macrostate, taking into account the energy levels and respective degenerate states of all the molecular species that intervene in the reaction. The population for each energy level is obtained by relating it to the partition function as per Equations 77 to 79.

$$N_j^{\mathrm{A}}=N^{\mathrm{A}}{\displaystyle \frac{g_j^{\mathrm{A}}{e}^{-{\epsilon }_j^{\mathrm{A}}/kT}}{Q^{\mathrm{A}}}}$$

(77)

$$N_j^B=N^B{\displaystyle \frac{g_j^B{e}^{-{\epsilon }_j^B/kT}}{Q^B}}$$

(78)

$$N_j^{\mathrm{AB}}=N^{\mathrm{AB}}{\displaystyle \frac{g_j^{\mathrm{AB}}{e}^{-{\epsilon }_j^{\mathrm{AB}}/kT}}{Q^{\mathrm{AB}}}}$$

(79)

9.3.1. Equilibrium Constant

The law of mass action is obtained by relating the several species of reactants and products, based on the zero-point energy as shown in Equation 80. Nonetheless, it is not practical to measure the number of molecules of the system, so Equation 81 relates the partial fraction of each species to the partition functions. The latter are directly dependent on volume, so the partial pressure fraction becomes solely a function of temperature, introducing the concept of a thermally variable equilibrium constant.

$${\displaystyle \frac{N^{\mathrm{A}}{N}^{\mathrm{B}}}{N^{\mathrm{AB}}}}=e^{-\Delta {\epsilon }_o/kT}{\displaystyle \frac{Q^{\mathrm{A}}{Q}^{\mathrm{B}}}{Q^{\mathrm{AB}}}}$$

(80)

$${\displaystyle \frac{p_{\mathrm{A}}{p}_{\mathrm{B}}}{p_{\mathrm{AB}}}}={\displaystyle \frac{kT}{V}}{e}^{-\Delta {\epsilon }_o/kT}{\displaystyle \frac{Q^{\mathrm{A}}{Q}^{\mathrm{B}}}{Q^{\mathrm{AB}}}}=K_p\left(T\right)$$

(81)

The equilibrium constant is more precisely defined for a general chemical reaction in Equation 82, where $v_i$ represents the stoichiometric coefficient of the species i, with positive values for products and negative values for reactants.

$$K_p\left(T\right)\equiv \prod _i{p}_i^{{\nu }_i}$$

(82)

9.3.2. Equilibrium Calculation

Qualitatively speaking, the equilibrium composition depends on temperature and pressure. To implement a proper algorithm capable of solving the equilibrium of a gaseous mixture, all relevant chemical reactions and their respective equilibrium constants must be defined. As a rule of thumb, for a mixture with x species containing y different elements, x-y independent chemical reactions must be considered. It’s likely to have more molecular species than equations, so other concepts need to be taken into account, such as conservation of atomic nuclei or Dalton’s law for gas partial pressure.

9.3.3. Equilibrium Enthalpy

Equation 83 represents the absolute enthalpy of a given species i, relating sensible enthalpy to the zero-point energy. Figure 13 illustrates the concept. The total absolute enthalpy of the entire mixture can be given by Equation 84, taking into account the mole-mass ratio ${\eta }_i$.

$$H_i={(H-E_o)}_i+E_{o_i}$$

(83)

$$h=\sum {\eta }_i{H}_i$$

(84)

Unfortunately, absolute enthalpy cannot be measured, a fact that does not undermine chemical studies, as the pertinent concept is the variation of enthalpy. In a chemically reacting mixture, it is necessary to establish a reference level for enthalpy, relative to which all energies are measured.

Equation 85 establishes the enthalpy variation for a system, with the second term relating to the heat of formation. This heat represents the energy released or absorbed during the formation of a given species from its constituent elements (for example, the natural element for species containing oxygen is $O_2$) under standard conditions ($T_s=298K$).

$$\Delta h=\Delta h_{\mathrm{sens}}+\Delta e_o$$

(85)

The heat of standard formation can be defined at absolute zero, where the temperature of both reactants and products is $0K$. In chemical reactions, the difference between the zero-point energy of the products and the reactants equals the difference between the heats of formation at $T=0K$. The heat of formation at absolute zero is considered the "effective" zero-point energy.

Equation 86 defines the total enthalpy, with the first term representing the sensible enthalpy calculated from statistical mechanics, and the second term reflecting the effective zero-point energy derived from the heat of formation at absolute zero, which has been experimentally determined and tabulated.

$$h=\sum _i{\eta }_i{(H-E_o)}_i+\sum _i{\eta }_i{(\Delta H_f)}_i^o$$

(86)

9.4. Nonequilibrium Systems

Chemical processes take time to occur and depend on molecular collisions, meaning processes are not instantaneous. The molecular collision frequency Z, proportional to $p/\sqrt{T}$, represents the number of collisions per second per molecule.

For example, an $O_2$ molecule needs around 20,000 collisions before becoming vibrationally excited. More broadly, the number of collisions depends on the molecule and its kinetic energy (i.e., its temperature). Increasing the temperature of the system leads to more violent collisions, which may provoke dissociation. For example, the $O_2$ molecule needs around 200,000 collisions to dissociate into atomic oxygen.

A system in chemical nonequilibrium simply hasn’t had enough time for all the necessary molecular collisions to occur and achieve a constant composition, pressure, and temperature at its chemical equilibrium. For instance, right downstream from a shock, the flow will have a non-equilibrium region due to pressure and temperature changes.

9.4.1. Vibrational Rate Equations

The expected number of transitions of the vibrational mode from level i to level $i+1$ is defined as the transition probability $P_{i,i+1}<<1$. Its product with the molecular collision frequency Z gives the number of transitions per particle per second, represented by $k_{i+1,i}$. Further multiplying by the number of molecules in the system gives the number of transitions from vibrational level i to $i+1$ per second.

Often, a given energy level will receive and give molecules more likely to its neighboring energy levels. So, the population of the energy level ${ϵ}_i$ will vary accordingly with Equation 87.

$${\displaystyle \frac{d{N}_i}{dt}}=k_{i+1,i}{N}_{i+1}+k_{i-1,i}{N}_{i-1}-k_{i,i+1}{N}_i-k_{i,i-1}{N}_i$$

(87)

Studying nonequilibrium conditions deepens the understanding of chemical equilibrium. For instance, the principle of detailed balancing states that in chemical equilibrium, the number of transitions from ${ϵ}_i$ to ${ϵ}_{i+1}$ and vice versa needs to be equal, in fact $d{N}_i/dt=0$.

In a nonequilibrium situation, an equilibrium vibrational energy can be established so the system moves in that direction. However, as the vibrational energy decreases, it is transferred to translational energy, increasing the temperature, which, in turn, raises the expected equilibrium vibrational energy to which the system is converging.

Equation 88 allows for the modeling of the vibrational energy as the system moves in the direction of chemical equilibrium. It’s important to note that, besides the current vibrational energy $e_{vib}$, the expected equilibrium vibrational energy $e_{vib}^{eq}$ is also a function of time. $\tau$ represents the model’s relaxation factor, which depends on the system’s temperature and pressure and combines the transition probability $P_{i,i+1}$ and collision frequency Z.

$${\displaystyle \frac{d{e}_{\mathrm{vib}}}{dt}}={\displaystyle \frac{1}{\tau }}(e_{\mathrm{vib}}^{\mathrm{eq}}-e_{\mathrm{vib}})$$

(88)

9.4.2. Rate of Chemical Reactions

A given chemical reaction will have two equilibrium constants, $k_f$ and $k_b$, referring to the forward and backward directions of the chemical equation. The rate of change in concentration of the reactant participant $X_i$ for the forward and backward reactions is given by Equations 89 and 90, respectively. Equation 91 represents the net rate in the system.

$${\displaystyle \frac{d\left[X_i\right]}{dt}}=(v_i^{\prime \prime }-v_i^{\prime })k_f\prod _i{\left[X_i\right]}^{{\nu }_i^{\prime }}$$

(89)

$${\displaystyle \frac{d|X_i|}{dt}}=-({\nu }_i^{\prime \prime }-{\nu }_i^{\prime })k_b\prod _i{\left[X_i\right]}^{v_i^{\prime \prime }}$$

(90)

$${\displaystyle \frac{d|X_i|}{dt}}=({\nu }_i^{\prime \prime }-{\nu }_i^{\prime })\left\{k_f\prod _i{\left[X_i\right]}^{{\nu }_i^{\prime }}-k_b\prod _i{\left[X_i\right]}^{{\nu }_i^{\prime \prime }}\right\}$$

(91)

The entire chemical equation can be modeled with an equilibrium constant K defined in Equation 92, with the prefix c if referring to species concentration or p if referring to species partial pressure. The equilibrium constants for most generic chemical equations have already been experimentally found.

$$K_c={\displaystyle \frac{1}{\mathcal{R}T}}{K}_p={\displaystyle \frac{k_f}{k_b}}$$

(92)

9.5. Modeling of High Enthalpy Flow

Chemical equilibrium denotes a constant chemical composition of the system, with consistent partial pressure of gaseous species. Thermodynamic equilibrium occurs when the system lacks pressure, temperature, velocity, and species concentration gradients, indicating that the energy levels are populated in accordance with the Boltzmann distribution.

When a system is in both chemical and thermodynamic equilibrium, it is referred to as being in complete thermodynamic equilibrium. However, in practical terms, this state is only achieved in a stationary gas and does not accurately represent flows.

It’s crucial to note that the Mach number loses its relevance in high enthalpy flows and, in nonequilibrium models, it forfeits its physical significance. Additionally, specific heats exhibit significant oscillations in value, making them less commonly used.

9.5.1. Local Equilibrium

The flow expanding in a rocket nozzle typically does not achieve complete thermodynamic equilibrium. However, if pressure and temperature gradients are sufficiently small, infinitesimal equilibrium can be assumed.

In this context, equilibrium is calculated for each point in the flow based on local pressure and temperature conditions, implying an infinite chemical rate of change.

For instance, through a shock, the temperature increase can lead to the gas becoming vibrationally activated and chemically reactive. To calculate the flow properties after the shock, the following iterative method can be employed:

• Assume an initial value for ${\rho }_1/{\rho }_2$ ( $0.1$ is often used);
• Calculate $P_2$ with Equation 93;
• Calculate $h_2$ with Equation 94;
• Calculate ${\rho }_2$ as a function of $P_2$ and $h_2$ using equations of state;
• Calculate the new ${\rho }_1/{\rho }_2$ ;
• Return to step 2. In case $P_2$ and $h_2$ have converge, follow to step 7;
• With $P_2$ and $h_2$ calculate $T_2$ using equations of state;
• Calculate $v_2$ with Equation 95.

$$p_2=p_1+{\rho }_1{u}_1^2\left(1-{\displaystyle \frac{{\rho }_1}{{\rho }_2}}\right)$$

(93)

$$h_2=h_1+{\displaystyle \frac{u_1^2}{2}}\left[1-{\left({\displaystyle \frac{{\rho }_1}{{\rho }_2}}\right)}^2\right]$$

(94)

$$u_2={\displaystyle \frac{{\rho }_1{u}_1}{{\rho }_2}}$$

(95)

In a calorically perfect gas, the conditions downstream of the shock depend solely on the upstream Mach number $M_1$. However, in a chemically reacting gas, properties depend on the local equilibrium established by $P_2$ and $T_2$, and consequently on the upstream temperature $T_1$, pressure $P_1$, and velocity magnitude $u_1$. For a thermally perfect gas, the pressure dependence is disregarded.

Equation 96 presents the equilibrium velocity of sound, which is established under the assumption that pressure and temperature within the shock wave remain constant and in local equilibrium.

$$a_e^2=\gamma RT{\displaystyle \frac{[1+(1/p){(\partial e/\partial v)}_T]}{[1-\rho {(\partial h/\partial p)}_T]}}$$

(96)

9.5.2. Frozen Flow

If chemical rates are considered null, the gas composition remains constant, allowing for analogies to a calorically perfect gas. However, this method can only be applied under the assumption of inviscid flow, as diffusion could alter the local gas composition.

A vibrationally frozen flow can be established, resulting in lower predicted temperatures since frozen energy is not transferred to other energetic modes. If the flow is also chemically frozen, specific heats remain constant. Equations 97 and 98 represent the specific heats as the sum of the frozen specific heat and the contribution due to chemical reactions. The additional contribution from chemical reactions is often dominant in terms of magnitude.

$$c_p=\sum _i{c}_i{c}_{p_i}+\sum _i{h}_i{\left({\displaystyle \frac{\partial c_i}{\partial T}}\right)}_p=c_{p_f}+\sum _i{h}_i{\left({\displaystyle \frac{\partial c_i}{\partial T}}\right)}_p$$

(97)

$$c_v=c_{v_f}+\sum _i{e}_i{\left({\displaystyle \frac{\partial c_i}{\partial T}}\right)}_v$$

(98)

Equivalent properties can be employed, such as considering the velocity at which sound would travel if the flow were not chemically reacting (frozen speed of sound, $a_f$) or the adiabatic index to plot an isentropic expansion (effective adiabatic index ${\gamma }_{eff}$). Although this approach improves results, they are still deficient in accuracy compared to numerical solutions.

9.5.3. Nonequilibrium Numerical Modulation

So, for chemical reactions with a finite, non-zero rate, numerical solutions to the governing equations must be employed. Chemical composition becomes a function not only of pressure, temperature, and velocity magnitude but also of the chemical reaction rate and duct geometry.

The continuity equation is split for each chemical species in the flow, as shown in Equation 99. Here, ${\dot{w}}_i$ represents the local rate of change of the species i due to chemical reactions, reflecting Equation 91. If viscosity is considered, a diffusive term must be added. The other governing equations undergo similar considerations.

$${\displaystyle \frac{\partial {\rho }_i}{\partial t}}+\nabla ·\left({\rho }_i\mathbf{V}\right)={\dot{w}}_i$$

(99)

9.5.4. High Enthalpy Shock

In a shockwave, the conversion of the flow’s kinetic energy into molecular internal energy occurs. In a calorically perfect gas, all of this energy is delivered to translational and rotational energy modes, resulting in a higher predicted temperature, $T_2$.

However, in a thermally perfect gas, the internal energy is distributed among not only translational and rotational modes but also vibrational modes, and even zero-point energies in the case of chemically reacting gas. This means that $T_2$ will be lower compared to the calorically perfect gas case.

Pressure, on the other hand, is less sensitive to these energy distributions as it is a "mechanical" property, depending mostly on the flow conditions themselves rather than on thermodynamics. A higher pressure upstream, $P_1$, will result in a higher pressure downstream, $P_2$, which can retard phenomena such as dissociation and ionization, leading to a higher downstream temperature as translational and rotational modes of energy receive more energy.

The shock front itself is so thin that it can always be considered frozen, and nonequilibrium conditions begin immediately downstream of the shock front.

9.5.5. High Enthalpy Quasi-1D Flow

In a nozzle, where the flow is likely to be chemically reactive, a quasi-1D analysis can be applied using local equilibrium, as changes in pressure and temperature are quickly corrected by the flow. This analysis resembles that of a calorically perfect gas, but with an iterative method to calculate the equilibrium. Pressure and temperature across the nozzle are functions of stagnation conditions and local velocity magnitude.

Figure 14 shows that in the reservoir, high temperatures lead to molecular dissociation, with atoms recombining through the expansion, as temperature decreases. In the expansion, temperature decreases less then with calorically perfect gas, as molecular recombination gives away energy due to zero-point energy decreasing.

If nonequilibrium is being considered, the sonic line, from local equilibrium considerations, will be slightly curved and downstream of the geometric throat [6].

10. Transonic Solution

The throat region of a convergent-divergent nozzle must be carefully studied, as it is where the flow chokes, and most contour design methods derive crucial boundary conditions. Early on, it was discovered that the sonic line is not straight but actually slightly convex in the upstream direction, a consequence of the nature intrinsic to the local flow, which exhibits both subsonic and supersonic phenomena, characterized by pockets at different regimes. It became of utmost importance to develop a sufficiently accurate method to model the velocity vectors of the flow in the throat region.

Hall [8] proposed the first method of a transonic solution near a de Laval nozzle, that is still applied by some contemporary literature. It consists of an expansion in order of $R_c^{-1}$, where $R_c$ is the radius of the contour immediately downstream of the geometric throat. By expanding the series of velocity to $R_c^{-3}$, parabolic, hyperbolic, and circular arc throats are considered. This method can be implemented for other divergent shapes with some adaptations.

Hall’s method demonstrates considerable accuracy compared to experimental results for $\gamma =1.4$, $R_c=5$ and two-dimensional flow. The series expansion to $R_c^{-1}$ yielded an error of just $0.5\%,$ and when fully expanded to $R_c^{-3}$, the error further decreased to only $0.02\%$.

For $R_c<1.5$, Hall’s solution starts to detriorate, and for values $R_c<1$, category of most of the current rocket nozzles, the method diverges in an oscillatory manner. Some inaccuracies can also be detected in the terms of the third degree parcel of the expansion for $u_3$ and $v_3$. To address these problems, Kliegel [9] presented a similar method, but the series are expansions of ${(R_c+1)}^{-1}$. This upgraded method can even consider sharp corner corresponding to $R_c=0$.

Hall’s solution employs cylindrical coordinates, but Kliegel defends that toroidal coordinates must be used as that’s the only way to exactly solve the wall boundary condition. Kliegel’s method shows good correlation with experimental results for $R_c=0.625$.

Nonetheless, Kliegel’s method still presents some errors. Specifically, in bringing the problem back to cylindrical coordinates from toroidal ones, it does not ensure that the equations of motion are respected in the cylindrical coordinates.

Dutton [10] addresses this problem with a ${(R_c+\eta )}^{-1}$ expansion, where $\eta$ is arbitrary. Dutton method respects the differential equations of motion for cylindrical coordinates when bringing the control volume back from toroidal. In fact, as seen in Figure 15 a), Dutton’s method is the most accurate for $R_c=0.625$ compared to experimental results.

An algorithm with Dutton’s method was implemented, presented in Appendix B with the respective Equations. Its sonic line for $R_c=0.625$ was traced and can be seen in Figure 15 b), where’s evident the in-house code provides a curve with similar shape and intercepts the axis and contour in the same points as Dutton’s graphic.

11. Method of Characteristics

The quasi-1D approach fails to consider the multidimensional of velocity and other flow properties, making it unfit to provide contour designs and analyses. Within a flow, the full-velocity potential equation takes the form of a partial differential equation (PDE). The Method of Characteristics (MOC) is a mathematical tool to solve PDEs by transforming them into ordinary differential equations (ODE).

To accomplish this, the MOC employs characteristic lines, inherent to the mathematical definition of a PDE, along which the PDE simplifies into an ODE. These lines maintain a mathematical link, facilitating the solving of the flow at different points, starting from known conditions, often derived from assumptions specific to a given nozzle design.

To apply the following description of the MOC, the flow must be assumed to be supersonic, steady, inviscid (with neglected boundary layer effects) and irrotational[6].

The potential velocity in Equation 100 is obtained by combining the continuity equation and Euler’s equation for a two-dimensional, irrotational flow. This results in a nonlinear PDE that becomes a hyperbolic equation if Equation 101 is satisfied, a condition met in supersonic flow.

$$\left[1-{\displaystyle \frac{{\Phi }_x^2}{a^2}}\right]{\Phi }_{xx}+\left[1-{\displaystyle \frac{{\Phi }_y^2}{a^2}}\right]{\Phi }_{yy}-{\displaystyle \frac{2{\Phi }_x{\Phi }_y}{a^2}}{\Phi }_{xy}=0$$

(100)

$${\displaystyle \frac{{\Phi }_x^2+{\Phi }_y^2}{a^2}}>1$$

(101)

The direction lines defining the hyperbola are known as characteristics, and in the context of the velocity potential in supersonic flows, coincide with the Mach lines.

As a function of both x and y, the velocity potential’s derivatives are expressed by Equations 102 and 103.

$$d\left({\Phi }_x\right)=du={\Phi }_{xx}dx+{\Phi }_{xy}dy$$

(102)

$$d\left({\Phi }_y\right)=dv={\Phi }_{xy}dx+{\Phi }_{yy}dy$$

(103)

By combining Equations 100, 102 and 103 in a system of linear algebraic equations with variables ${\Phi }_{xx}$, ${\Phi }_{yy}$ and ${\Phi }_{xy}$, and solving for ${\Phi }_{xy}$, using Crammer’s, Equation 104 can be obtained.

$${\Phi }_{xy}={\displaystyle \frac{\left|\begin{array}{ccc}1-\frac{u^2}{a^2}& 0& 1-\frac{v^2}{a^2}\\ dx& du& 0\\ 0& dv& dy\end{array}\right|}{\left|\begin{array}{ccc}1-\frac{u^2}{a^2}& \frac{-2uv}{a}& 1-\frac{v^2}{a^2}\\ dx& dy& 0\\ 0& dx& dy\end{array}\right|}}={\displaystyle \frac{Numerator}{Denominator}}$$

(104)

At a given point in the flow, ${\Phi }_{xy}$ implies that along a direction $dy/dx$, the velocity undergoes changes corresponding to the respective $du$ and $dv$ components, except for a specific direction, where the denominator equals zero. Along this direction, the derivative lacks a specific finite value, posing an inconsistency. Therefore, for this scenario to be indeterminate, the numerator must also be zero.

From this mathematical definition, despite the flow’s properties being continuous in the absence of shocks, the first derivative must be indeterminate along a characteristic line. This condition allows for the crossing of several streamlines.

Setting the denominator of Equation 104 to zero, results in Equation 105. The slope of the characteristic lines that pass through a given point can then be obtained using Equation 106, simplified into Equation 107 [11].

$$\left[1-{\displaystyle \frac{u^2}{a^2}}\right]{\left({\displaystyle \frac{dy}{dx}}\right)}_{char}^2+\left[1-{\displaystyle \frac{v^2}{a^2}}\right]+{\displaystyle \frac{2uv}{a^2}}{\left({\displaystyle \frac{dy}{dx}}\right)}_{char}^2=0$$

(105)

$${\left({\displaystyle \frac{dy}{dx}}\right)}_{char}={\displaystyle \frac{\frac{-uv}{a^2}\pm \sqrt{\frac{u^2+v^2}{a^2}-1}}{1-\frac{u^2}{a^2}}}$$

(106)

$${\left({\displaystyle \frac{dy}{dx}}\right)}_{char}=tan(\theta \pm \mu )$$

(107)

Considering Equation 107, it becomes evident that, in supersonic flow, any given point is intersected by two characteristic lines. In the case of sonic flow, a single characteristic line exists, and its slope equals the Mach angle. For subsonic flow, the Mach angle becomes imaginary, signifying that the characteristic lines take an elliptical form [6].

Examining Figure 16, two characteristics are observed passing through a designated point A: a left characteristic $C_{+}$ with a slope of $tan(\theta +\mu )$ and a right characteristic $C_{-}$ with a slope of $tan(\theta -\mu )$

To solve the flow along a characteristic line, the governing PDE describing the flow reduces to ODEs, also known as compatibility equations. These equations are obtained by setting the determinant of the numerator in Equation 104 to zero, resulting in Equation 108.

$${\displaystyle \frac{dv}{du}}=-{\displaystyle \frac{1-\frac{u^2}{a^2}}{1-\frac{v^2}{a^2}}}{\left({\displaystyle \frac{dy}{dx}}\right)}_{char}$$

(108)

Equation 109 represents the set of ODEs obtained from Equation 108. With algebraic manipulation it leads to the expressions illustrated in Equations 110 and 111 for a planar flow.

$$d\theta =\pm \sqrt{M^2-1}{\displaystyle \frac{dV}{V}}$$

(109)

$$\theta +\nu \left(M\right)=c^{te.}=K_{-}along{C}_{-}$$

(110)

$$\theta -\nu \left(M\right)=c^{te.}=K_{+}along{C}_{+}$$

(111)

In cases where the flow is unknown but the characteristic constants have been established at the nodes, Equations 112 and 113 can be employed to determine the flow axial deviation and Prandtl-Meyer angle, respectively [11].

$$\theta ={\displaystyle \frac{1}{2}}\left[\left(K_{-}\right)+\left(K_{+}\right)\right]$$

(112)

$$\nu \left(M\right)={\displaystyle \frac{1}{2}}\left[\left(K_{-}\right)-\left(K_{+}\right)\right]$$

(113)

Rocket nozzles commonly exhibit axisymmetric behavior, making cylindrical coordinates a more appropriate choice for describing the flow. The prior description of the MOC can be modified, with the key distinction lying in the compatibility equations denoted in Equations 114 and 115. In the axisymmetric MOC, the assumption of a constant relation along a given characteristic line no longer holds. Consequently, the relation between the variables $\theta$ and $\nu$ requires integration, which can be achieved by employing the finite differences method described in Appendix C accordingly to Zebbiche [12] and Restrepo [13].

$$d(\theta +\nu )={\displaystyle \frac{1}{\sqrt{M^2-1}-cot\theta }}{\displaystyle \frac{dr}{r}}={\displaystyle \frac{sin\mu sin\theta }{sin(\theta -\mu )}}{\displaystyle \frac{dr}{r}}along{C}_{-}$$

(114)

$$d(\theta -\nu )={\displaystyle \frac{1}{\sqrt{M^2-1}+cot\theta }}{\displaystyle \frac{dr}{r}}=-{\displaystyle \frac{sin\mu sin\theta }{sin(\theta +\mu )}}{\displaystyle \frac{dr}{r}}along{C}_{+}$$

(115)

For a planar context, initiating from established boundary conditions, characteristic lines radiate to form a grid. The MOC is initially employed to solve the flow at each node. Subsequently, the slopes obtained are utilized to construct the grid, imparting geometric dimensions to the nodes, solving the flow. To note that, when approximating a characteristic line between two nodes i and i+1, to a straight line, the sloop should be an average as depicted in Equation 116. For an axisymmetric case, both the nodes and their coordinates are solved at once and one node at the time as there’s no longer a constant through a characteristic line.

$$tan\left[{\displaystyle \frac{1}{2}}({\theta }_i+{\theta }_{i+1})\pm {\displaystyle \frac{1}{2}}({\mu }_i+{\mu }_{i+1})\right]$$

(116)

If the MOC is employed to draw the contour of a nozzle, there are two main philosophies: either drawing a line that follows the orientation of the flow at each point or setting mass conservation to draw a specific streamline. The origin of the grid depends on the particular nozzle design and associated assumptions.

A great programming language to employed MOC algorithms is Python, as it is one of the most widely used in aerospace engineering, including by NASA, especially in non-safety critical systems [14]. It is characterized as a high-level programming language that emphasizes readability and encourages the creation of open-source libraries. Additionally, Python is free to use and does not require a compilation step, which accelerates the edit-test-debug cycle and prevents segmentation faults from incorrect user usage through the use of exceptions.

12. CFD - Computational Fluids Dynamics

Computational Fluid Dynamics (CFD) is a branch of fluid mechanics, including any numerical method that can analyze and solve problems that involve flows. By discretizing the fluid domain into a grid and solving governing equations using numerical methods, CFD enables the visualization and understanding of complex fluid behaviors.

This powerful tool is extensively used in engineering and physics, allowing for the prediction of aerodynamics, heat transfer, and other fluid-related phenomena. CFD simulations play a vital role in optimizing designs, evaluating performance, and reducing the reliance on costly physical prototypes in various industries.

Although the MOC can be considered part of CFD, this section of the book aims to provide a broader understanding of CFD, incorporating simplifications typical of supersonic flows.

Reference [6] serves as the basis for this Section, except the last Subsection.

12.1. Differential Conservation Equations for lnviscid Flows

Traditional, fluid analysis, and consequently conservation equations, are set for a control volume $\mathcal{V}$ delimited by surface area S. So, for any vector function A and scalar function $\Phi$Equations 117 and 118 are true for any control volume, respectively.

$${∯}_{S}A·dS=∰(\nabla ·A)d\mathcal{V}$$

(117)

$${∯}_S\Phi dS={∰}_{\mathcal{V}}(\nabla \Phi )d\mathcal{V}$$

(118)

12.1.1. Continuity

Equation 119 represents the continuity equation for a given control volume. From Equation 117 with $A=\rho V$, Equation 120 is obtained.

$$-{∯}_S\rho V·dS={∰}_{\mathcal{V}}{\displaystyle \frac{\partial \rho }{\partial t}}d\mathcal{V}$$

(119)

$${∰}_{\mathcal{V}}\left[{\displaystyle \frac{\partial \rho }{\partial t}}+\nabla ·\left(\rho V\right)\right]d\mathcal{V}=0$$

(120)

For Equation 120 to be null, the first parcel of the integral could be annulled by the second. This is impractical, as the analysis must be true for any random control volume, so the only chance for it to be zero is if the control volume is a single point, being precisely the goal in CFD. So Equation 121 is obtained, representing the differential equation of continuity.

$${\displaystyle \frac{\partial \rho }{\partial t}}+\nabla ·\left(\rho V\right)=0$$

(121)

12.1.2. Momentum

Equation 122 represents the momentum conservation for a given control volume. If combined with Equation 118 with $\Phi =p$, Equation 123 is obtained. Equation 124 represents its scalar decomposition in the x direction.

$${∰}_{\mathcal{V}}\rho fd\mathcal{V}-{∯}_{S}pdS={∰}_{\mathcal{V}}{\displaystyle \frac{(\partial \rho V)}{\partial t}}d\mathcal{V}+{∯}_S(\rho V·dS)V$$

(122)

$${∰}_{\mathcal{V}}\rho fd\mathcal{V}-{∰}_{\mathcal{V}}\nabla pd\mathcal{V}={∰}_{\mathcal{V}}{\displaystyle \frac{(\partial \rho V)}{\partial t}}d\mathcal{V}+{∯}_S(\rho V·dS)V$$

(123)

$${∰}_{\mathcal{V}}\rho f_{x}d\mathcal{V}-{∰}_{\mathcal{V}}{\displaystyle \frac{\partial p}{\partial x}}d\mathcal{V}={∰}_{\mathcal{V}}{\displaystyle \frac{(\partial \rho u)}{\partial t}}d\mathcal{V}+{∯}_S(\rho V·dS)u$$

(124)

If Equation 117 assumes $A=\rho V$, Equation 125 is obtained. Further consideration of Equation 124 leads to a single integral equal to zero, as depicted in Equation 126.

$${∯}_S(\rho V·dS)u={∯}_S\left(\rho uV\right)·dS={∰}_{\mathcal{V}}\nabla ·\left(\rho uV\right)d\mathcal{V}$$

(125)

$${∰}_{\mathcal{V}}\left[\rho f_x-{\displaystyle \frac{\partial p}{\partial x}}-{\displaystyle \frac{\partial \left(\rho u\right)}{\partial t}}-\nabla ·\left(\rho uV\right)\right]d\mathcal{V}=0$$

(126)

Solving Equation 126 with the same philosophy of the continuity, Equations 127 , 128, 129 represent differential equations of momentum conservation for each Cartesian direction.

$${\displaystyle \frac{\partial \left(\rho u\right)}{\partial t}}+\nabla ·\left(\rho uV\right)=-{\displaystyle \frac{\partial p}{\partial x}}+\rho f_x$$

(127)

$${\displaystyle \frac{\partial \left(\rho v\right)}{\partial t}}+\nabla ·\left(\rho vV\right)=-{\displaystyle \frac{\partial p}{\partial y}}+\rho f_y$$

(128)

$${\displaystyle \frac{\partial \left(\rho w\right)}{\partial t}}+\nabla ·\left(\rho wV\right)=-{\displaystyle \frac{\partial p}{\partial z}}+\rho f_z$$

(129)

12.1.3. Energy

Equation 130 represents the energy equation for a given control volume. Combined with Equation 117 and setting $A=\rho (e+V^2/2)V$ and $A=eV$, Equation 131 is reached.

$$\begin{array}{cc}\hfill & {∰}_{\mathcal{V}}\rho \dot{q}d\mathcal{V}-{∯}_{S}pV·dS+{∰}_{\mathcal{V}}\rho (f·V)d\mathcal{V}=\hfill \\ & {∰}_{\mathcal{V}}{\displaystyle \frac{\partial }{\partial t}}\left[\rho \left(e+{\displaystyle \frac{V^2}{2}}\right)\right]d\mathcal{V}+{∯}_S\rho \left(e+{\displaystyle \frac{V^2}{2}}\right)V·dS\hfill \end{array}$$

(130)

$${∰}_{\mathcal{V}}\left\{\rho \dot{q}-\nabla ·\left(pV\right)+\rho (f·V)-{\displaystyle \frac{\partial }{\partial t}}\left[\rho \left(e+{\displaystyle \frac{V^2}{2}}\right)\right]-\nabla ·\left[\rho \left(e+{\displaystyle \frac{V^2}{2}}\right)V\right]\right\}d\mathcal{V}=0$$

(131)

Solving Equation 131 with the previous philosophy, Equation 132 is set and represents the differential equations of energy conservation.

$${\displaystyle \frac{\partial }{\partial t}}\left[\rho \left(e+{\displaystyle \frac{V^2}{2}}\right)\right]+\nabla ·\left[\rho \left(e+{\displaystyle \frac{V^2}{2}}\right)V\right]=-\nabla ·\left(pV\right)+\rho \dot{q}+\rho (f·V)$$

(132)

12.2. Introduction to Finite Differences

In finite differences methods, flow’s properties are known and calculated in discrete nodes. This unit process can be better understood by looking at Figure 17. This method of solving the flow can be advantageous compared to the MOC, as sometimes characteristic lines can distort and following them is no longer possible, besides the impossibility of ensuring an uniform grid.

At each node, properties need to allow the expansion of modeling equations in the form of a Taylor’s series, as in Equation 133 for velocity $u_{i,j}$ along the abscissa axis. Homologous equations are derivated for other properties and Cartesian directions.

$$u_{i+1,j}=u_{i,j}+{\left({\displaystyle \frac{\partial u}{\partial x}}\right)}_{i,j}\Delta x+{\left({\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}\right)}_{i,j}{\displaystyle \frac{{(\Delta x)}^2}{2}}+\cdots$$

(133)

12.2.1. Types of Differences

When building a CFD model, the way ${\left({\displaystyle \frac{\partial u}{\partial x}}\right)}_{i,j}$ is found must be choose in accordance to the problem that’s being solved.

The forward difference gives the derivative as the linear gradient between the node and the immediate next downstream node as per Equation 134.

$${\left({\displaystyle \frac{\partial u}{\partial x}}\right)}_{i,j}={\displaystyle \frac{u_{i+1,j}-u_{i,j}}{\Delta x}}$$

(134)

The rearward difference gives the derivative as the linear gradient between the node and the immediate prior upstream node as per Equation 135. For this case, the Taylor’s expansion is given by Equation 136.

$${\left({\displaystyle \frac{\partial u}{\partial x}}\right)}_{i,j}={\displaystyle \frac{u_{i,j}-u_{i-1,j}}{\Delta x}}$$

(135)

$$u_{i+1,j}=u_{i,j}+{\left({\displaystyle \frac{\partial u}{\partial x}}\right)}_{i,j}(-\Delta x)+{\left({\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}\right)}_{i,j}{\displaystyle \frac{{(-\Delta x)}^2}{2}}+\cdots$$

(136)

The central difference gives the derivative as the linear gradient between the prior upstream node and the immediate next downstream node as per Equation 137.

$${\left({\displaystyle \frac{\partial u}{\partial x}}\right)}_{i,j}={\displaystyle \frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}}$$

(137)

This way, the partial derivatives in the flow’s ruling equations can be replaced by algebraic differences.

12.2.2. Order of Accuracy

The truncation of the Taylor’s series introduces some errors, that will increment with the increasing of the partitions $\Delta x$ and/or $\Delta y$. The more terms that are used, the more accurate the solution will be, but also lot more of computational power will be required. Factors like convergence behavior and stability must also be considered.

Usually first-order accuracy is employed with the Taylor serie’s given by Equation 138.

$$u_{i+1,j}=u_{i,j}+{\left({\displaystyle \frac{\partial u}{\partial x}}\right)}_{i,j}\Delta x+\cdots$$

(138)

If a more accurate solution is desired, second-order accuracy can be employed as suggested in Equations 133 or 136. For viscous flows, second-order derivatives are present in the momentum and energy equations.

The following discussion will focus on first-order accuracy for inviscid supersonic flow.

12.2.3. Difference Equations

When the partial differentials in the flow’s equations are replaced with the differences presented before, difference flow’s equations are obtained.

For example, Equation 121 for steady 2D flow becomes Equation 139. Letting $F=\rho u$ and $G=\rho v$ and applying a forward difference in x and a central difference in y, Equation 140 is obtained, representing the difference equation.

$${\displaystyle \frac{\partial \left(\rho u\right)}{\partial x}}+{\displaystyle \frac{\partial \left(\rho v\right)}{\partial y}}=0$$

(139)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{F_{i+1,j}-F_{i,j}}{\Delta x}}={\displaystyle \frac{G_{i,j+1}-G_{i,j-1}}{2\Delta y}}\hfill \\ & \iff F_{i+1,j}=F_{i,j}+{\displaystyle \frac{\Delta x}{2\Delta y}}(G_{i,j+1}-G_{i,j-1})\hfill \end{array}$$

(140)

The governing equations can be expressed in a generic form for 3D steady flow, as shown in Equation 141. In this equation, each vector is defined, with the first term reflecting the continuity equation, the last term representing the energy equation, and the middle three terms each corresponding to a Cartesian projection of the momentum equation.

$${\displaystyle \frac{\partial F}{\partial x}}+{\displaystyle \frac{\partial G}{\partial y}}+{\displaystyle \frac{\partial H}{\partial z}}+J=0$$

(141)

$$F=\left\{\begin{array}{c}\rho u\\ \rho u^2+p\\ \rho vu\\ \rho wu\\ \rho (e+V^2/2)u+pu\end{array}\right\}$$

(142)

$$G=\left\{\begin{array}{c}\rho v\\ \rho uv\\ \rho v^2+p\\ \rho wu\\ \rho (e+V^2/2)v+pv\end{array}\right\}$$

(143)

$$H=\left\{\begin{array}{c}\rho w\\ \rho uw\\ \rho vw\\ \rho w^2+p\\ \rho (e+V^2/2)w+pw\end{array}\right\}$$

(144)

$$J=\left\{\begin{array}{c}0\\ \rho f_x\\ \rho f_y\\ \rho f_z\\ \rho \dot{q}+\rho (u{f}_x+v{f}_y+w{f}_z)\end{array}\right\}$$

(145)

12.2.4. Explicit v.s. Implicit Methods

Considering the grid of Figure 17, with the values of the vertical line i known. The downstream line i+1 values can be directly calculated from the values of line i.

This illustrates an explicit solver, where downstream nodes are calculated solely based on upstream information. Its the easiest to implement but the values of $\Delta x$ and $\Delta y$ must be carefully chosen so the solution is stable. It leads to more dense grids with higher computational cost.

Another algorithm may consider the variation of F as the average of the value of G of the upstream line i of known values and the unknown downstream line i+1, as suggested by Equation 146.

$$\begin{array}{cc}\hfill & {\displaystyle \frac{F_{i+1,j}-F_{i,j}}{\Delta x}}={\displaystyle \frac{1}{2}}\left[\left({\displaystyle \frac{G_{i,j+1}-G_{i,j-1}}{2\Delta y}}+{\displaystyle \frac{G_{i+1,j+1}-G_{i+1,j-1}}{2\Delta y}}\right)\right]\hfill \\ & \iff F_{i+1,j}=F_{i,j}+{\displaystyle \frac{\Delta x}{4\Delta y}}(G_{i,j+1}+G_{i+1,j+1}-G_{i,j-1}-G_{i+1,j-1})\hfill \end{array}$$

(146)

This way, a system of equations is built, so the solution does not only depend on upstream known flow but also on the downstream line that’s being calculated. It must be applied to several points at the same time so enough equations are available to build a proper solving matrix.

This implicit method has the advantage of permitting a more stable solutions for more scarce grids, allowing for faster advancement of the algorithm, with less computational costs, although the solving matrix takes some of the computational gains. Nonetheless, this type of method is harder to implement.

12.3. MacCormack’s Technique

There are several finite differences schemes with the present Subsection focusing on the one proposed by Robert MacCormack at the NASA Ames Research Center, that was very popular in the 70’s and 80’s. Although it currently has been surpassed by modern schemes, it is considered student friendly and is worth exploring to have a concrete idea of what this kind of algorithms do.

MacCormack’s method is applied to supersonic inviscid steady flow that has hyperbolic equations. It has two steps giving it a second-order accuracy although Taylor series are truncated at the first degree.

Starting by considering Equation 141 for 2D flow with no source in the fashion of Equation 147, the partial derivative can be given as an average between the known upstream line i and the unknown downstream line i+1 as suggested in Equation 148

$${\displaystyle \frac{\partial F}{\partial x}}=-{\displaystyle \frac{\partial G}{\partial y}}$$

(147)

$$F_{i+1,j}=F_{i,j}+{\left({\displaystyle \frac{\partial F}{\partial x}}\right)}_{\mathrm{ave}}\Delta x$$

(148)

Firstly, a predictor step applies a forward difference to calculate the predicted ${\overline{F}}_{i+1,j}$ value as in Equation 149 for all the points in the upstream line i.

$$\begin{array}{cc}\hfill & {\overline{F}}_{i+1,j}=F_{i,j}+{\left({\displaystyle \frac{\partial F}{\partial x}}\right)}_{i,j}\Delta x\hfill \\ & with{\left({\displaystyle \frac{\partial F}{\partial x}}\right)}_{i,j}=-{\left({\displaystyle \frac{\partial G}{\partial y}}\right)}_{i,j}=-{\displaystyle \frac{G_{i,j+1}-G_{i,j}}{\Delta r}}\hfill \end{array}$$

(149)

Then, in a corrector step , and after ${\overline{F}}_{i+1,j}$ is used to calculate ${\overline{G}}_{i+1,j}$, a reward difference is used in Equation 150 to find the average partial differential that is then calculated in Equation 151

$${\left({\displaystyle \frac{\overline{\partial F}}{\partial x}}\right)}_{i+1,j}=-{\displaystyle \frac{{\overline{G}}_{i+1,j}-{\overline{G}}_{i+1,j-1}}{\Delta y}}$$

(150)

$${\left({\displaystyle \frac{\partial F}{\partial x}}\right)}_{ave}={\displaystyle \frac{1}{2}}\left[{\left({\displaystyle \frac{\partial F}{\partial x}}\right)}_{i,j}+{\left({\displaystyle \frac{\overline{\partial F}}{\partial x}}\right)}_{i+1,j}\right]$$

(151)

Homologous thinking can be applied to solve the flow in another Cartesian directions.

There are still some applications of the MacCormack’s technique in contemporary literature like Li [15] that proved a good correlation with experimental results, inclining in modeling an oblique shock.

12.4. Boundary Conditions

Figure 18 highlights some challenges of finite differences, namely flow’s boundary conditions. Abbett’s method consist in an algorithm for inviscid, steady and supersonic flow that applies MacCormack’s technique but only with forward differences for both prediction-correction for physical walls, like in point 2.

The previous process will results in calculated velocity, $v_{calc}$, that is not tangent to the surface. So a correction must be set as in Equation 152 that is reflected in all flow properties, meaning an isentropic expansion or compression takes place, deviating the flow an angle $\theta$.

$$v_{act}=v_{calc}+\theta$$

(152)

Abbett’s method can also be applied for nodes downstream of shocks, like point 4, with the prediction-correction consisting on solely forward or reward differences. After the oblique shock, the actual angle of the flow is easily given by Equation 22.

12.5. Stability Criterion

As it has been refereed before, for explicit solvers, $\Delta x/\Delta y$ must be choose so the algorithm is stable in predicting the solution. Larger values of the previous ratio leads to excessive truncation errors decreasing accuracy.

For instances, it can be argued that, for a vertical line i, the node $i+1,j$ should be upstream of the interception of the right and left characteristics of nodes $i,j+1$ and $i,j-1$, respectively. Courant-Friedrichs-Lewy (CFL) criterion reflects precisely this and is mathematically formalized in Equation 153.

$$\Delta x\le {\displaystyle \frac{\Delta y}{{|tan(\theta \pm \mu )|}_{max}}}$$

(153)

12.6. Shock Modeling

A shock-fitting approach, seen in Figure 19 a) can be implement when the shock localization and angle are known, modeling it as discontinuity in flow’s properties.

A shock-capturing approach, seen in Figure 19 b) will find the shock, with the grid covering all the flow and taking free stream properties as boundaries. The shock will be the region with strong properties gradients. Although there’s no prior need of knowing the shock, the grid needs to be denser and some calculations at the far-field are yield that will not be relevant for the final solution

12.7. Turbulence Models

The occurrence of laminar or turbulent flow is dependent on the relative influence of viscous and inertial effects, as described by the Reynolds number. As the Reynolds number increases, the development of a laminar flow into a turbulent flow becomes more likely.

Turbulent flow is characterized by the existence of eddies that create large fluctuations in velocity and other properties, as depicted in Figure 20. Reynolds’ decomposition states that the value of a velocity component at a certain point is the sum of its average value and a fluctuation that is a function of time, as seen in Equations 154 and 155. It is important to note that turbulent motion is random, and therefore a statistical treatment is essential [16].

$$u\left(t\right)=\overline{u}+u^{\prime }\left(t\right)$$

(154)

$$v\left(t\right)=\overline{v}+v^{\prime }\left(t\right)$$

(155)

In summary, turbulent flow is irregular, highly diffusive and dissipative, three-dimensional, and has a large Reynolds number. The eddies that characterize turbulence transfer their kinetic energy to smaller eddies in a process called the energy cascade. This continues until the eddies become so small that the viscous stresses are enough to dissipate the kinetic energy into internal energy. These small scales are also called Kolmogorov scales and are much smaller than the characteristic length of the flow. They can be approximated with isotropic properties, although the flow is still considered a continuum medium.

To model turbulence in CFD, Reynolds decomposition must be applied to the governing equations, such as the Reynolds-Averaged Navier-Stokes (RANS) equations. In RANS and other governing equations, the values of properties such as velocity, pressure and other scalars are replaced by their average, leaving a term called the Reynolds tensor, $\overline{u_i^{\prime }{u}_j^{\prime }}$, with the fluctuations. Turbulence models propose ways to handle the Reynolds tensor by adding more equations and variables to the physical model, reflecting assumptions [17].

The Reynolds tensor can be related to the mean velocity gradients through Equation 156 by adding a new variable, ${\mu }_t$, called turbulent viscosity, which is dependent on the flow characteristics rather than the fluid conditions, unlike dynamic viscosity. This modulation of the Reynolds tensor is known as the Boussinesq Approach, and the turbulence models presented in this Subsection are introduced in order to compute the turbulent viscosity. It offers a relatively low computational cost since it assumes the turbulent viscosity to be an isotropic scalar, being sufficiently accurate for wall boundary layers, mixing layers, jets, etc.., where the flow’s shear is dominated by only one of the turbulent shear stresses.

$$-\rho \overline{u_i^{\prime }{u}_j^{\prime }}={\mu }_t\left({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}\right)-{\displaystyle \frac{2}{3}}\left(\rho k+{\mu }_t{\displaystyle \frac{\partial u_k}{\partial x_k}}\right){\delta }_{ij}$$

(156)

The Spalart-Allmaras model only adds one transport equation, representing the turbulent viscosity. It was specially developed for aerospace applications, namely, wall-bounded flows, with accurate results for adverse pressure gradients. Nonetheless, it shows larger errors for some free shear flows, especially plane and round jet flows.

The $k-ϵ$ model is a high-Reynolds number model and adds two equations: one to model the turbulence kinetic energy, k, and another to model the turbulence dissipation rate, $ϵ$. Very favored in industrial flow and heat transfer simulations, this semi-empirical model provides a fair accuracy for a wide range of flows with an economic usage of computational resources. In its standard form, it assumes that the flow is fully turbulent, ignoring the effects of molecular viscosity.

The RNG $k-ϵ$ model is more accurate and reliable for a wider class of flows, and derives from the standard version by adding a statistical technique called renormalization group theory (RNG). Some of the improvements are explained by: an extra term in the turbulence dissipation rate equation capable of improving the accuracy for rapidly strained flows, considering swirl effects, instead of user-specified and constant turbulent Prandtl numbers the statistics gives an analytical formula for it, and an analytically-derived differential formula for effective viscosity providing means to simulate low-Reynolds number flows.

The realizable $k-ϵ$ model is a relatively recent approach that satisfies mathematical constraints associated with Reynolds stresses. Therefore, this model offers improved correspondence to the physical phenomena observed in turbulent flows. It achieves this by employing an alternative modulation of turbulent viscosity and dissipation rate. The dissipation rate is obtained from an exact equation for the transport of mean-square vorticity fluctuations. Similar to the RNG $k-ϵ$ model, the realizable $k-ϵ$ model demonstrates superior performance, particularly in flows characterized by strong streamline curvature, vortices and rotation.

Lastly, the $k-\omega$ model is an empirical model and adds two equations: one to model the turbulence kinetic energy, k, and another to model the specific turbulence dissipation rate, $\omega$. It suits low-Reynodls numbers, compressibility and shear flow spreading. Originally it was freestream sensitive but has been improved over the years so it can predict free shear flows [18].

12.8. Turbulence Models for Nozzle Simulation

Typical simulations of the flowfield inside a nozzle focus on properly modeling the flow and boundary layer near the contour and the jet’s interaction with ambient air. Turbulence models are employed to close the differential form of governing equations. Since high Reynolds numbers are expected, large eddy simulation fails, and RANS equations are used, necessitating the development of turbulence models to predict scalar transport terms and Reynolds stresses.

Due to high property gradients, meshes are dense, making the computational cost of second-order RANS prohibitive. First-order RANS turbulent models like mixing lengths and Spalart-Allmaras are inadequate for capturing turbulent effects in detail within a nozzle’s flowfield. Two-equation models, based on Boussinesq’s concept of viscous turbulence ${\mu }_t$, are commonly adopted. Most models include equations for turbulent kinetic energy (k) and either k dissipation ($ϵ$) or k specific dissipation rate ($\omega$).

One widely used turbulence model for nozzle flow is the shear stress transport (SST) k-$\omega$, which calculates turbulent viscosity considering the transport effects of the main shear stress. Its greatest advantage lies in gradually transitioning within the boundary layer from a traditional k-$\omega$ to a version of k-$ϵ$, with the latter having lower computational cost and being applied to more mesh cells in practical terms.

The SST k-$\omega$ model has proven effective for ducted flows from injectors, accurately predicting shock waves and boundary layer detachment, albeit with slight overprediction of shear stresses, potentially leading to higher velocities in secondary flow. Detailed mathematical modeling of SST k-$\omega$ can be found in Reference [19].

13. Modern Methods for Nozzle Design

The method of characteristics proves to be a reliable tool for designing supersonic nozzles, mainly due to its simplicity. However, it comes at the cost of neglecting some important factors relevant in nozzles for wind tunnels and/or high enthalpy flow. For example, by considering an inviscid flow, the Method of Characteristics neglects real gas effects like changing specific heats and the boundary layer.

Considering local vibrational equilibrium or applying a displacement to the contour representing the boundary layer are simply implemented methods to improve the Method of Characteristics. However, they still lack important features specially for hypersonic flow, characterized by high enthalpy effects and where a parallel uniform exit flow is mandatory.

Current computational power allows for more complex methods to design a supersonic nozzle, mostly by providing relatively fast tools that modulate the flow and permit subtle changes to the contour. These methods account for more sophisticated approaches to real flow.

13.1. Improved Method of Characteristics - Changing Adiabatic Index

More accurate nozzle contours can be obtained by considering real gas effects, which can be easily integrated into the MOC algorithm. In high enthalpy flows, real gas effects become evident, and the calorically perfect assumption induces significant errors.

An algorithm of the MOC can be implemented, but with the local vibrational equilibrium assumption with a thermically perfect gas, where the adiabatic index is now a function of temperature as per Equation 157, with $\Theta =3055.556K$ [20].

$$\gamma \left(T\right)={\displaystyle \frac{{\gamma }_{perf}-1}{1+({\gamma }_{perf}-1)\left[{\left(\frac{\Theta }{T}\right)}^2\frac{e^{\Theta /T}}{{(e^{\Theta /T}-1)}^2}\right]}}$$

(157)

For rocket nozzles this method can provide a more reliable contour, as so for supersonic wind tunnels [21]. However, for hypersonic wind-tunnels other approaches are needed, as for $M>8$ dissociation becomes relevant, with the thermically perfect gas delivering inaccurate results.

13.2. Improved Method of Characteristics - Boundary Layer Correction

While the perfect gas assumption simplifies calculations, it neglects the crucial influence of the boundary layer. This thin region near the nozzle wall, where viscous forces are significant, impacts the flow behavior. Fluid in contact with the wall adheres due to the no-slip condition, resulting in zero velocity there. As the fluid moves farther from the wall, its velocity gradually increases, reaching the undisturbed freestream velocity outside the boundary layer.

The boundary layer’s thickness is influenced by the fluid’s viscosity, object shape, and freestream velocity. Generally, thinner layers occur for low-viscosity fluids and streamlined objects [22].

The presence of a boundary layer introduces a displacement thickness, as shown in Equation 158. This effectively reduces the nozzle’s cross-sectional area, leading to lower exit Mach numbers compared to those predicted solely by the geometrical area ratio.

$${\delta }^{*}={\int }_0^{\delta }\left(1-{\displaystyle \frac{{\rho }_e{u}_e}{{\rho }_{BL}{u}_{BL}}}\right)$$

(158)

Furthermore, the flow field and boundary layer exhibit a two-way interaction, impacting the final exit velocity profile. Higher velocities and stagnation temperatures lead to a thicker boundary layer, making its influence non-negligible in high-enthalpy flows. Consequently, the inviscid contour requires displacement, as detailed in Equation 159[23].

$$\left\{\begin{array}{cc}\hfill x_{BLcorrected}& =x_{\mathrm{inviscid}}-{\delta }^{*}sin\theta \hfill \\ \hfill y_{BLcorrected}& =y_{\mathrm{inviscid}}+{\delta }^{*}cos\theta \hfill \end{array}\right.$$

(159)

There are simple algorithms available to predict the boundary layer, and CFD can be utilized for iterating the final contour. Two essential parameters characterizing boundary layers are the shape factor H, which indicates whether the boundary layer is expected to be laminar or turbulent, and the momentum thickness $\phi$, defined in relation to the momentum flow rate within the boundary layer. These three parameters are interrelated as described by Equation 160.

$$H={\displaystyle \frac{{\delta }^{*}}{\phi }}$$

(160)

The accuracy of modeling a boundary layer relies on understanding the velocity profile within the boundary layer. Equation 161 presents the von Kármán momentum integral equation for planar compressible flow.

$${\displaystyle \frac{C_f}{2}}={\displaystyle \frac{\mathrm{d}\phi }{\mathrm{d}x}}+{\displaystyle \frac{\theta }{U_{\infty }}}{\displaystyle \frac{\mathrm{d}{U}_{\infty }}{\mathrm{d}x}}(2+{\displaystyle \frac{{\delta }^{*}}{\phi }})={\displaystyle \frac{d\phi }{dx}}+\phi \left[{\displaystyle \frac{2-M_e+H}{M_e(1+\frac{\gamma -1}{2}{M}_e^2)}}{\displaystyle \frac{d{M}_e}{dx}}\right]$$

(161)

Various simplified methods can be employed to solve the momentum differential equation. For instance, applying a Stewartson transformation to H and $\phi$ and substituting them into Equation 161. The resulting equation can then be solved using a fourth-order Runge-Kutta method. This approach involves calculating several flow variables, including the friction coefficient $C_f$ and incompressible equivalents.

Another approach involves using turbulent boundary-layer theory for a flat plate, which can be a reasonable approximation for a 2D nozzle, despite the Mach number not being constant. The momentum equation is simplified to Equation 162, and the boundary-layer displacement thickness is determined by Equation 163. The Reynolds number is defined by Equation 164, and viscosity is calculated using Sutherland’s formula [24] for nitrogen gas, as shown in Equation 165.

$${\displaystyle \frac{d\phi }{dx}}={\displaystyle \frac{C_f}{2}}$$

(162)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{{\delta }^{*}}{x}}\approx 0.046{\mathrm{Re}}_x^{-1/5},{\mathrm{Re}}_x<{10}^7\hfill \\ \hfill & {\displaystyle \frac{{\delta }^{*}}{x}}\approx 0.018{\mathrm{Re}}_x^{-1/7},{10}^5<{\mathrm{Re}}_x<{10}^9\hfill \end{array}$$

(163)

$$Re={\displaystyle \frac{\rho VD}{\mu }}={\displaystyle \frac{P_e{M}_{e}y}{{\mu }_e}}\sqrt{{\displaystyle \frac{\gamma }{R{T}_e}}}$$

(164)

$${\displaystyle \frac{\mu }{{\mu }_0}}={\left({\displaystyle \frac{T}{T_0}}\right)}^{3/2}{\displaystyle \frac{T_0+S_{\mu }}{T+S_{\mu }}}\iff {\mu }_{e_{N_2}}=1.663\times {10}^{-5}{\left({\displaystyle \frac{T_e}{273}}\right)}^{3/2}{\displaystyle \frac{273+107}{T_e+107}}$$

(165)

Ding [25] employed both methods to simulate the boundary-layer displacement of an engine inlet of a hypersonic vehicle with axisymmetric geometry. The study found good accuracy in the numerical integration of the momentum differential equation with CFD experiments. However, the simplification of considering turbulent boundary layers over flat plates under-predicts the displacement by around $35\%$.

13.2.1. Limitations of Method of Characteristics based Contours

For high enthalpy nozzles, particularly hypersonic nozzles operating at Mach numbers higher than 8, the previous methods are not precise enough. On one hand, the flow becomes chemically reactive, and on the other hand, the boundary-layer thickness becomes on the order of the size of the exit radius. Irregular flow rates arising from compressive waves at inflection points could be addressed by implementing Parabolic Navier-Stokes equations in the design method.

Furthermore, when nozzles are longer, the boundary-layer and main flow may interact, so the characteristic line is not reflected at the frontier between the two regions but within the boundary layer. This invalidates the previous method, as characteristic lines are deviated, not completely canceling the expansion waves [26].

13.2.2. CFD based Methods

The MOC may provide an accurate initial approach when designing a supersonic nozzle, but it is very limited in taking into account real gas effects.

CFD has been developed over the years with various solvers that allow for the consideration of real gas effects, from boundary-layer wall modeling to chemically reacting gases.

Design methods are proposed, often starting with a contour developed through the MOC, and iteratively obtaining the flowfield of the nozzle and changing the contour to approach design parameters such as the exit Mach number or exit flow angle.

However, conducting a flow simulation in CFD for each iteration, along with the need for sophisticated optimization algorithms, makes these methods computationally expensive.

For example, [27] designs a minimum-length nozzle using an Euler solver for each contour, employing a surrogate-based optimization algorithm to find the nozzle contour with the optimal thrust coefficient. The results showed little difference compared to a reduced computational cost method that employs a free-form deformation parameterization of the contour, opening doors for the development of shape optimization strategies, specially for preliminary design.

14. Method of Characteristics for Conventional Nozzles Design

15. Conical Nozzle

As seen in Figure 21, a conical nozzle is the simplest types of divergent nozzle contour. It is characterized by the throat diameter $R_t$, exit diameter $R_e$, and divergence angle $\alpha$. Due to its practical manufacturing process, has historically been one of the most used contours [2]. The bigger the divergence angle, the shorter the nozzle will be for a given area ratio, as this is the design desire, since length is a fair indicator of a nozzle’s weight.

In theory, the flow velocity at the exit could be considered one-dimensional, covering all the exit’s area, but due to the wall having an angle right before the exit, not all velocity vectors are parallel to the rocket, reducing the effective thrust and can even generate side loads and[2]. The divergence loss of traction can be estimated with Equation 166.

$$\lambda ={\displaystyle \frac{1+cos\alpha }{2}}$$

(166)

This type of nozzle is very susceptible to flow separation and overexpansion at lower altitudes. The thrust losses increase with the decrease of nozzle length by manipulation of the divergence angle. Taking into account nozzle’s length, performance and weight, an half-angle of 15º offers the best compromise. Literature states it should never surpassed 25º degrees as separation becomes to intense for practical applications.

Nonetheless, when a low area ration is desired, and fabrication simplicity is privileged in deterioration of aerodynamic performance, which are design parameters for some missiles, there’s no need for more complex designs and conical nozzles offer a good trade-off [28].

The emergence of CFD methods for nozzle analysis has left a gap in the past decades regarding the application of the MOC to solve conical nozzles. Despite the lack of recent research, a MOC algorithm was implemented to gain a deeper understanding of certain phenomena associated with conical nozzles. The algorithm solves the flow inside a prescribed two-dimensional expansion ramp contour, defined by the half-angle, the flow’s adiabatic index, and the desired Mach number at the exit, assuming quasi-1D expansion.

Besides the pre-defined conditions for applying the MOC, from its mathematical formulation, it was also assumed that the flow inside a conical nozzle behaves in accordance with Figure 22, where a centered expansion fan at the lip "a" of the throat is responsible for expanding the flow from $\theta =0^{\circ }$ to $\theta =\alpha$. If the flow at point "c" has $\theta =\alpha$ and $M=M_{exit}$, its left characteristic with $K_{+}=\alpha -\nu \left(M_{exit}\right)$ can be traced until it meets the axis at point "b," where $K_{+}=0-{\nu }_b$, leading to ${\nu }_b=\nu \left(M_{exit}\right)$. From point "b," the right characteristic $K_{-}=0+\nu \left(M_{exit}\right)$ is traced back to the centered expansion fan at point "a," where, due to the nature of the referred structure, ${\nu }_a={\theta }_a$. Writing the equation of the last right characteristic at "a," which is the same for point "b," it can be concluded that the maximum expansion angle, which must equal the cone half-angle, should be given by Equation 167.

$${\theta }_{max}=\alpha ={\displaystyle \frac{\nu \left(M_{exit}\right)}{3}}$$

(167)

Unfortunately, the obtained results prove the algorithm fails to capture the assumed behavior. Firstly, the conditions expected for point "c" were observed inside the nozzle, as depicted in Figure 23. This discrepancy could be attributed to the divergence angle that the flow will experience at the exit. Thus, two quasi-1D expansions were considered from points "c" and "b". For the specific case of $M_{exit}=2.4$ and $\gamma =7/5$, the average value between the Mach number at the exit lip and center-line was equal to the design Mach number. Regrettably, this observation did not hold for other design conditions, with the averaged exit Mach number being higher than the design value and increasing with the deviation from the previous design conditions.

If the ramp half angle is smaller than the one calculated in Equation 167, the expansion fan will fail to provide the necessary acceleration. Consequently, secondary expansion waves must be present along the contour wall. However, this can lead to isentropic losses because the flow should not curve more than the half angle of the ramp.

In the case of a ramp half angle "$\alpha$ larger than the one calculated in Equation 167, the initial expansion may be excessive, and the characteristic lines from it seem to influence the flow outside the nozzle. At the exit plane, shocks may occur in order to raise the pressure of a most likely overexpanded flow.

In reality, conical nozzles lack an isentropic solution, even when operating at design conditions. This phenomenon is associated with two factors. First, due to the expansion occurring along the cone, the flow expands inwards, resulting in the inevitable appearance of a shock at the exit plane. Secondly, even for small half-cone angles, the throat is too sharp, always causing an oblique shock inside the nozzle. This creates an observable pattern of a double Mach diamond, with a thickness peaking near design operation and a small Mach disc inside the nozzle. This is visualized in Figure 24 where it must be highlighted that even at design operation, a Mach disk is found inside the nozzle.

Shadowgraph images reveal that underexpanding conditions cause an increase in expansion at the throat, potentially making the throat shock to overtake the Mach diamond structure. Increasing the exit Mach number results in more elongated cells at the Mach diamond, with the two shocks converging and coalescing after the initial pair of compression and expansion waves.

In cases of overexpansion, decreases in Mach number also lead to the coalescence of the two Mach diamonds, primarily due to the formation of a Mach disk caused by the lip shock, with a diameter close to half of the exit [29]. Similar to other traditional nozzles, in overexpansion, the entire contour is not fully utilized, leading to separation at the end. This scenario poses risks due to the generation of side loads, especially during start-up [30] and reattachment during the transition to total contour utilization.

Historically, to mitigate the presence of the throat shock and achieve an isentropic conical nozzle, the suggestion of incorporating an initial expanding arc has been made. Yet, utilizing the asymmetrical MOC reveals that a weak oblique shock will always consistently emanate from the junction between the arc and the cone unless the junction is optimized to render the shock and its effects almost imperceptible [31].

In conclusion, conical nozzles never operate in the third critical point, always incurring isentropic losses in addition to divergence. Furthermore, some engineering challenges are posed such as side loads, which can compromise high reliability and mission optimization.

16. Bell Nozzle

As the divergence angle of a conical nozzle is increased, the flow divergence at the exit can increase significantly. To minimize this effect, more complex nozzle designs are necessary. For example, the bell nozzle is a more complex design that can provide a better performance.

Figure 25 shows a comparison between an equivalent conical nozzle and two bell nozzles with different lengths and contour angles. The bell nozzles are shorter, with lengths of 80% and 60% of the length of the conical nozzle, and have smaller contour angles at the exit, with values of ${11}^{\circ }$ and 8.$5^{\circ }$ compared to the 15º of the conical nozzle. These design changes help to maintain a more axial flow at the exit, reducing the degree of divergence and improving the overall performance of the nozzle.

The gradual reverse of the divergence angle through the nozzle is possible by a ${20}^{\circ }{-50}^{\circ }$ divergence right after the nozzle, which does not introduce significant losses because of the huge favorable pressure gradient that permits this rapid expansion without separation.

This makes the bell nozzle the most widely used nozzle design in modern rocketry and from which most advanced designs are derived, as it already shows good performance and reliability results [2].

The flow structure inside a bell nozzle can be modeled with the MOC, considering an initial kernel region immediately after the throat, where most of the expansion occurs. In this region, a thick grid is necessary to capture the property gradients. Then, a straightening section follows, where the characteristics are mostly constant from the kernel to the contour, serving as the means of constructing the latest.

16.1. Minimum Length Contour

16.1.1. Minimum Length Contour - Two-Dimensional

The most straightforward design for a bell nozzle with a zero-degree inclination at the exit takes the assumption that all expansion takes place at a centered expansion fan located at the throat’s lip. The role of the contour is to align the resultant flow without introducing isentropic losses.

The preceding assumption asserts that the left characteristic, uniting nodes 34 and 35 in Figure 26, is defined by a constant $\theta =0$ and $\nu =\nu \left(M_e\right)$. Node 34 represents the intersection of the last right characteristic with $K_{-}=\nu \left(M_e\right)$, originating from the lip at "a," where, owing to the centered expansion fan, $\theta =\nu$. Utilizing Equation 110, this results in ${\theta }_a={\theta }_{MAX}=0.5\nu \left(M_e\right)$.

Therefore, to generate the Minimum Length contour, right characteristic lines are traced from the lip "a", in accordance with the expanding flow, ensuring that $0<\theta \le {\theta }_{MAX}$, until the nozzle’s axis is intersected and undergo reflection. Subsequently, the points on the last right characteristic line are projected until crossing the contour. This contour is a straight line with a slope equal to the average flow inclination at the previous point and the point being projected.

The algorithm described above has been implemented in an in-house Python code and validated using values from Reference [6]. The grids, illustrated in Figure 26 and Figure 27, are almost identical, with flow parameters at the nodes matching within at least two decimal points. The exception is the Mach angle, which may differ, likely due to the use of different methods for Prandtl-Meyer angle inversion. Notably, the in-house code exhibits a more accurate area ratio.

Mishra [32] also employs a similar algorithm, calculating the exit Mach number with the MOC resulting in 3.8759. Nevertheless, when the contour is tested in Ansys Fluent, the measured value is only 3.4, suggesting that the difference raises from the MOC’s omission of considerations for isentropic losses and the boundary layer.

The exploration of the Python code is presented in the Tables of Appendix H. As anticipated, the area ratio and contour length increase with higher exit Mach numbers and a decrease in the adiabatic index. The algorithm’s accuracy, when compared to quasi-1D flow, improves with smaller Mach numbers and higher adiabatic indexes. For higher Mach numbers and lower adiabatic indexes, the area ratio and length of the resulting 2D nozzle are more sensitive to the number of characteristic lines employed. It was observed that after considering 50 characteristic lines, the area ratio and nozzle length exhibit minimal change for any input pair. Further increases in the number of lines are only justifiable when a higher-resolution contour is required.

When compared to a straight 2D expansion ramp, with a constant 8º angle, it becomes evident that, for smaller Mach numbers, the contoured nozzle is actually larger than the expansion ramp. So, for shorter nozzles, the combination of an expansion fan followed by a straightening section to achieve total axial flow at the exit results in a longer nozzle, doubling the length in some cases. In reality, accepting a certain degree of divergence angle at the exit can induce minimal thrust losses, which are compensated for by a significantly shorter nozzle.

16.1.2. Minimum Length Contour - Two-Dimensional with Changing Adiabatic Index

The in-house algorithm was modified to consider $\gamma$ as a function of temperature. An iterative method was implemented, utilizing the Mach number of the respective node. The process assumes isentropic expansion from stagnation with $\gamma ={\gamma }_{iteration}$, providing a temperature that is used to calculate $\gamma$ for the next iteration. In the first iteration, ${\gamma }_0={\gamma }_{perf}=1.4$. Special attention needs to be given to the Prandtl-Meyer angle inversion, as it also becomes an iterative process.

The dimensions of some nozzles designed with the improved in-house code can be found in Appendix I.

Since it’s assumed a thermally perfect gas, only stagnation temperature was changed, with the algorithm having a constant stagnation pressure at the chamber.

Larger values of stagnation temperature led to longer nozzles, as a higher temperature represents lower $\gamma$ values, which from calorically perfect gas studies showed that higher area ratios and longer nozzles are needed. There are some exceptions, as there’s a shorter nozzle for $M_{\mathrm{exit}}=4.5$ at a stagnation temperature of $1500K$.

It must be highlighted that some assumptions keep holding true, like considering the centered wave expansion at the throat lip ends at an angle equal to half of the exit Prandtl-Meyer angle, as the measured exit Mach number is quite accurate.

For the same stagnation temperature, a bigger Mach number at the exit has a larger $\gamma$ at the exit, as static temperature is lower.

It must be concluded that even with longer projected nozzles, the contour is, in theory, more reliable. Even when the stagnation temperature is $500;K$ for $M_{\mathrm{exit}}=4.5$ with ${\gamma }_{\mathrm{exit}}=1.4$, the dimensions of the contour are slightly different from those designed with the perfect gas assumption, showing some influence in the kernel from the thermally perfect gas.

Sun [33] utilized JANNAF’s database to implement a frozen flow with varying specific heats for nozzle contour design, modifying Rao’s method. The modified contour, when tested in CFD, achieved up to $1\%$ additional specific impulse at the design point and demonstrated superior performance in off-design conditions.

16.1.3. Minimum Length Contour - Boundary-Layer Correction

Although the boundary-layer theory for a flat plane is less accurate, due to its simplicity, it was the method chosen, as it provides an understanding of the boundary layer, aligning with the goals of the present study.

Two in-house codes were developed, one assuming calorically perfect gas and the other assuming thermically perfect gas, to implement a boundary-layer correction to the inviscid contour. Results from varying the design Mach number and stagnation properties can be found in Appendix J. Figure 28 illustrates the inviscid and BL-corrected contours, both for thermically perfect gas, with parameters set to $M_{\mathrm{exit}}=2.4$, $T_t=1500;K$, $P_t=1;MPa$, and 200 characteristic lines.

The analysis revealed that, for all cases, the thermically perfect gas assumption resulted in a thicker boundary layer due to the lower adiabatic indexes, allowing for lower Reynolds numbers, as per Equation 164.

An increase in the stagnation temperature of the main flow was found to increase the boundary-layer thickness, while decreasing the stagnation pressure produced a similar effect.

The application of the boundary layer displacement thickness method models the development of the boundary layer, and ${\delta }^{*}$ at the exit is larger for longer nozzles, even with a higher Reynolds number. A simplification of this method involves neglecting the boundary layer at the throat, which can be substantial and crucial for designing a nozzle with the proper area ratio.

16.1.4. Minimum Length Contour - Axisymmetry

It is also possible to construct the contour of an axisymmetric nozzle with the previous minimum length assumption adapted to the new flow’s geometry. There is no longer a constant along a characteristic line, so the flow is solved in the fashion described in Appendix C.

Characteristic lines still emanate from lip A in accordance with an expansion fan, but the maximum angle ${\theta }^{*}$ is not defined a priori. Characteristics are traced with increasing angles, and the algorithm checks if the desired exit Mach number is achieved every time a new characteristic intercepts the axis, allowing point B to be found.

Then, from point B, a left characteristic is projected with a fixed increasing x, assuming a constant Mach number ($M=M_E$) and angle ($\theta =0^{\circ }$), which will create a second grid between it and the last right characteristic from the kernel, by expanding the left characteristics.

Subsequently, the contour starts from lip A with a slope equal to ${\theta }^{*}$ until it meets the extension of the kernel’s left characteristics. From that point, the contour assumes the angle of the flow at that point of the characteristic, and the process repeats until reaching the left characteristic of point B, defining point E. It’s worth noting that the properties between two nodes, in this second grid, are assumed to change linearly with respect to the radius y, since a sufficient number of characteristic lines are employed.

An example of the characteristics grids for $M=2.4$, $\gamma =1.4$ and 7 lines can be seen on Figure as so the generated contour Figure 30.

When compared with the contours published in the literature, the in-house code yields a similar result, as seen in Figure 31. Some of the differences in the curve can be attributed to the interpolation of contour points and their calculation, especially in the approximation of property changes between nodes. Nonetheless, both contours have approximately the same length and area ratio.

The code was tested similarly to that for two-dimensional flow, and the results are found in Appendix K. The findings are comparable to those of two-dimensional flow. One difference arises from the fact that, in almost all cases, the minimum length bell nozzle is longer than a corresponding quasi-1D conical nozzle, making parabolic and truncated nozzles more interesting subjects of study. This trend decreases, tho, for higher Mach numbers at the exit or lower adiabatic indexes, where in both cases there’s a need for more expansion.

Due to the way the algorithm finds ${\theta }^{*}$, when a low number of characteristics is considered, the exit Mach number is higher than the design one, varying directly with the found ${\theta }^{*}$. This suggestion is backed by the exit Mach number and the value of ${\theta }^{*}$ itself. As gases with lower $\gamma$ have less aggressive expansions, the results tend to be more accurate.

When a considerable number of characteristics are taken into account, the algorithm always predicts an area ratio lower than the theoretical quasi-1D model, which can be attributed to the way the contour is traced, by employing linear relations to solve some unknown nodes and projecting the contour’s streamline. There’s a trend for the length and area-ratio to stabilize, with the number of characteristics, and is closer to the quasi-1D analysis for high $\gamma$. The length of the nozzle increases if the exit Mach number calculated is higher than the Mach number assumed at design.

It can be concluded that the minimum length philosophy for bell nozzles, especially if axisymmetric, is not only more challenging, time-consuming and complex than conical nozzles but also often results in longer and heavier nozzles. Anyhow, the previous sentence needs to be carefully interpreted, as further studies involving CFD and even experimental research need to be conducted for the entire flight envelope to assess and potentially discard of these bell nozzles. Conical nozzles, never have an isentropic solution, adding losses in addition to divergence. Since having a totally axial flow at the exit is not acceptable for rocket applications, the developed bell contours can, for instance, be truncated to save weight, with only slight divergence at the exit causing minimal thrust loss.

16.2. Parabolic Contour

A parabolic nozzle consists of a bell-shaped contour where the initial expansion section is defined by a arc-circle, followed by a straightening section traced as a parabola, as suggested in Figure 32. The design parameters include the exit Mach number, the radius of the initial expansion arc-circle R, and the respective arc-circle’s angle, coincident with the flow ${\theta }_N$, and the exit divergence angle ${\theta }_E$. Additionally, the equivalent conical nozzle, and the respective fraction f can be inputs, as it will be explored latter.

Rao [34] proposed a method to obtain the most efficient nozzle contour for a prescribed length and area ratio. Essentially, the method involved finding the initial arc-circle that would promote optimum expansion by employing the MOC and determining the zero of the first derivative of an integral representing the thrust function, subject to certain constraints, using the Lagrangian multiplier method. The contour would then be traced by mass conservation along the characteristic lines in the straightening section. It was later discovered that this method can be simplified into tracing a parabola with minimal thrust losses, so parabolic nozzles are often considered to have a Thrust Optimized Parabolic (TOP) contours.

An in-house code was developed to draw a parabolic bell nozzle following the method presented in Reference [35]. It begins by defining the exit Mach number, determining the area ratio $ϵ$ accordingly with Equation 30, and specifying the reference conical nozzle of half-angle $\alpha$ and length fraction f. The length and height of the parabolic nozzle are then calculated using Equations 168 and 169, respectively.

$$L_N=f\left({\displaystyle \frac{\left(\sqrt{ϵ}-1\right)R_t}{tan\left(\alpha \right)}}\right)$$

(168)

$$R_e=\sqrt{ϵ}{R}_t$$

(169)

Once , ${\theta }_N$ and ${\theta }_E$ are prescribed, the parabola can be drawn. A Bézier quadratic curve is suggested, and to draw it, three points are advised, as per Figure 33. Point N has the coordinates of the end of the arc-circle, and point E corresponds to the height and length of the nozzle. Point Q is the intersection of the straight lines that pass through points N and E, with ${\theta }_N$ and ${\theta }_E$ as slopes, respectively. The coordinates of the discreet points of the parabola are calculated with Equations 170 and 171. Figure 34 shows an example of the results of the algorithm for parabolic bell nozzle contours.

$$x\left(t\right)={(1-t)}^2{N}_x+2(1-t)t{Q}_x+t^2{E}_{x}0\le t\le 1$$

(170)

$$y\left(t\right)={(1-t)}^2{N}_y+2(1-t)t{Q}_y+t^2{E}_{y}0\le t\le 1$$

(171)

Another in-house algorithm was developed to analyze the performance of the drawn parabolas. A study was performed, varying the number of points in the initial line (mesh independence study), the exit Mach number, the arc-circle angle ${\theta }_N$, the exit angle ${\theta }_E$ and cone half angle $\alpha$ (which were considered constant), the fraction of the equivalent conical nozzle f, and chamber stagnation conditions. For the sake of simplicity, the variables were changed one at a time, basing the standard case defined in Table 1.

The MOC is initiated at a vertical line, set at the abscissa of the sonic point of the axisymmetry axis. The initial points’ flow properties are given by Dutton’s transonic solution. The kernel is calculated, and then the nodes on the last right characteristic progressively make their left characteristic intercept the contour, as in Case D of the axisymmetric MOC, forming a secondary grid. Once a left characteristic crosses the exit plane, the previous one is used to calculate the nozzle performance.

Between each discreet pair of consecutive nodes, linear properties variations with y were considered. Thrust is given by Equation 172 and the mass flow by Equation 173, where $\varphi$ is the inclination of the straight line that unites the two nodes. For the in-house code, an intermediary node is calculated for each right and left characteristic that crosses the exit plane, with the respective abscissa. Then, fractional thrust and mass flow are calculated between each of the new nodes. It is important to note that chamber conditions are given, and local pressure and temperature are determined by an isentropic expansion with $R=8314;J/kg.K$. For thrust calculation, the ambient pressure is set to ensure that the nozzle is operating on its design point. Lastly, the specific impulse is given by Equation 2. The results from the study that was carried out can be found in Appendix L.

$$\begin{array}{cc}\hfill & F={\int }_i^{i+1}\left[(p-p_{amb})+\rho V^2{\displaystyle \frac{sin(\varphi -\theta )cos\theta }{sin\varphi }}\right]2\pi ydy=\hfill \\ \hfill & \left[(p_{avg}-p_{amb})+{\rho }_{avg}{V}_{avg}^2{\displaystyle \frac{sin(\varphi -{\theta }_{avg})cos{\theta }_{avg}}{sin\varphi }}\right]\pi (y_i^2-y_{i+1}^2)\hfill \end{array}$$

(172)

$$\dot{m}={\int }_i^{i+1}\rho V{\displaystyle \frac{sin(\varphi -\theta )}{sin\varphi }}2\pi ydy={\rho }_{avg}{V}_{avg}{\displaystyle \frac{sin(\varphi -{\theta }_{avg})}{sin\varphi }}\pi (y_i^2-y_{i+1}^2)$$

(173)

It must be highlighted that, upon examining the formed grid from the in-house code for all cases, some of the right characteristics reflecting from the parabola coalesce, indicating the presence of a weak oblique shock, as shown in Figure 35. Since it is a weak shock, theoretically exhibiting almost isentropic behavior, no more sophisticated analyses than the MOC will be employed.

An average of the exit node’s Mach number and flow angle $\theta$ is calculated to better understand the velocity profile and explain variations in the estimated specific impulse. Usually, the averaged $\theta$ is higher than the exit angle due to the compression of the parabola, generating a bigger slope than the wall where there’s a higher node density. The exit Mach number is higher at the axis, which is reasonable as the characteristic that emanates from the exit lip intercepts the axis upstream of the axis interception with the exit, and some expansion still occurs along it.

Initially, a mesh independence study was carried out for the standard case by varying the number of nodes at the initial line. Little changes were observed for the specific impulse with 50 or more initial nodes. Nonetheless, the Mach number at the exit decreases with more points, approaching the design Mach number. This may happen due to a better modulation of the expansion. But since the average Mach number and $\theta$ are practically constant, this issue was overlooked. Henceforth, 200 initial nodes were used, expecting better modulation of the flow.

Designing parabolas with increasing exit Mach numbers, showed that, in theory, there’s an improvement in performance, even surpassing that of the equivalent conical nozzle. However, this needs to be interpreted with caution, as this comparison is based on a quasi-1D analysis of the conical nozzle, where the exit flow is considered to have a constant $\theta =\alpha$. In the present study, the average $\theta$ decreases with increasing Mach number, being lower than the cone half angle $\alpha$. Also, there’s higher coalescence of characteristic lines, so the oblique shock may no longer be safe to assume isentropic. $M_{exit}=2.4$ provides the relatively lowest and most accurate averaged Mach number, being higher for other cases. There’s also a lower limit for the acceptable design Mach number, as some short nozzles don’t have parabolic solutions.

The initial circle arc angle ${\theta }_N$ also plays an important role as it defines the initial expansion. Higher ${\theta }_N$ leads to higher specific impulse and lower exit Mach numbers. In theory, the oblique shock gets more evident due to the aggressive initial expansion and contour discontinuity. This may be the reason with the Mach numbers decreases and the specific impulse is compensated by higher thrust due to positive pressure differences. It also provides a lower average $\theta$, also responsible for net thrust increases.

Comparing with other conical nozzles, by equally changing ${\theta }_e$ and $\alpha$, it is seen that allowing bigger divergences at the exit of the parabolic contour leads to a loss of specific impulse. Lower design angles have even lower averaged $\theta$ and higher Mach numbers.

When the parabola is designed with a shorter length than the comparable conical nozzle, there are more losses due to divergence, caused by a more aggressive expansion that doesn’t allow the straightening of the flow and a higher Mach number that promotes negative pressure differences. The oblique shock is stronger for longer nozzles, though, as being closer to the length of the conical nozzle leads to an almost straight parabola, inducing a bigger discontinuity in the contour.

The increasing of the stagnation temperature leads to a bigger specific impulse. In the ideal gas, the Mach number at each node doesn’t change, so as temperature increases, local velocity does too, leading to an equal increase in mass flow and thrust. As losses of specific impulse, compared to quasi-1D analysis, are temperature independent, losses are the same, in fraction.

The ideal gas made the performance of the nozzle and the MOC grid itself stagnation pressure independent, especially since in all analyses, it was allowed the nozzles to operate at their design point.

Real gas analysis show the major advantage of a TOP contoured nozzle is that it can have a length that is a fraction of a conical nozzle, with the more common fraction being $80\%$, as going under $70\%$ results in isentropic losses and minimal thrust losses compared to other bell nozzles that don’t allow for divergence at the exit. For example, an $80\%$ parabolic nozzle produces $99.8\%$ of the specific impulse of a full-length bell [35]. The Korean KSLV-II rocket uses a TOP contour in its third stage for vacuum operation. The discontinuity between the arc circle and the parabola induces a shock that causes flow detachment and creates a trapped vortex, which causes particularly dangerous side loads during start-up as the flow is adhering to the nozzle. When the separation bubble bursts, the recirculation caused by backflow increases the local pressure, moving the separation point on the contour upstream [36].

The previous mechanism, responsible for side loads during engine start-up, can be observed in Figure 36, often referred to as the transition between Free Shock Separation (FSS) and Restricted Shock Separation (RSS). The danger arises from the fact that the mechanism does not work in an axisymmetric fashion, especially the reattachment of the flow in the RSS mode. The detachment typically originates at the end of the arc-circle due to the discontinuity, not in contour or slope but curvature [37].

While investigating Truncated Ideal Compressed (TIC) contoured bell nozzles, Takahashi [38] linked the transition between FSS and RSS to the interaction between the separation point and the Mach disk. Furthermore, in TIC nozzles, there’s an abrupt fall of pressure after the initial expansion followed by a slight increase. This changes the wall pressure, moving the separation point and causing it to interact with the Mach disk, leading to RSS situations with considerable side-loads.

Stark [39] experimentally studied both TIC and TOP nozzles to further understand the phenomena mentioned earlier, especially why the flow was being redirected to the wall, generating a RSS situation. For nozzles operating at higher pressure ratios, discontinuities in wall pressure were found but attributed to nitrogen condensation, highlighting a key phenomenon in real flow. CFD and experimental studies showed that the Mach disk was either convex or concave, but the convexity did not explain the RSS mode for overexpanding nozzles at low pressure ratios and did not seem to affect the intensity of side-loads. In some experimental trials, the RSS mode was observable for a pressure ratio of 5. The mechanism behind the transition was the separation shock that disturbed the Mach disk, tilting it and causing the flow to be redirected to the contour wall, leading to reattachment of the flow.

Bakulu [40] further explores the mechanisms behind the FSS to RSS transition. Acoustic studies suggested that transonic resonance may be present near the axis, between the Mach disk and the exit, for a TIC nozzle, generating a standing pressure wave downstream of the separation shock, capable of moving the separation point. Although the exact physical mechanism could not be explained, for a TIC nozzle, this transonic resonance can be responsible for up to $20\%$ of the side loads generated during start-up due to the transition from FSS to RSS.

16.3. Truncated Ideal Compressed Contour

The TIC nozzle, often comparable and sometimes confused with a TOP nozzle, is derived from an ideal contour. This ideal contour can be traced in a similar fashion to a minimum length nozzle, but with a arc-circle throat and the contour is traced with mass conservation considerations. It’s important to note that the design Mach number of the ideal nozzle must always be greater than that of the TIC nozzle.

The ideal contour is then truncated at the abscissa with the ordinate equal to the desired exit area ratio for the TIC nozzle. A compression factor, as per Equation 174, is determinable, and it is linearly applied to all points P’ of the truncated contour, as shown in Equation 175.

$$C=1-{\displaystyle \frac{x_e^{\prime \prime }-x_a}{x_e^{\prime }-x_a}}$$

(174)

$$x_{p^{\prime \prime }}=(x_{p^{\prime }}-x_a)(1-C{\displaystyle \frac{x_{p^{\prime }}-x_a}{x_{e^{\prime }}-x_a}})+x_a$$

(175)

The process from the ideal nozzle, to truncted nozzle and the compressed nozzle is better understood with Figure 37

The compression introduces a discontinuity in the slope of the curve at point A, as the angle $\theta$ at the beginning of the contour increases. To address this, a solution is to slightly move the contour downstream and expand the arc circle until both angles meet, ensuring the continuity of the first derivative of the curve. This also changes slight the radius ratio of the nozzle.

However, this correction comes with consequences. Now, the initial expansion will be even more aggressive, with the possibility of characteristic lines coalescing into a shock. This shock will increase local pressure and even wall pressure, potentially enhancing the produced thrust.

Hoffman [41] utilized the MOC to compare the performance of TIC and TOP nozzles. The study revealed that TOP nozzles always outperformed TIC, albeit by a small margin ranging from $0.34\%$ to $0.04\%$. This indicates that TIC is a viable alternative when TOP nozzles are not suitable. It was also found that longer initial ideal nozzles, meaning those with a higher area ratio, usually don’t have shocks inside when truncated and compressed. This is because they are less compressed, and the truncation almost leaves them with the desired length.

An in-house code was developed to design a TIC nozzle with a arc-circle with $R=0.625$ and ${\theta }_N=30º$. The last node led to an ideal contour with $M_{\mathrm{ideal}}=4.69$. It was truncated for $M_{\mathrm{exit}}=3.5$ and later compressed to half the length of a conical nozzle with $\alpha =15º$. 200 nodes were used in the initial line, similar to the code to evaluate the TOP contour performance. These three nozzles can be seen in Figure 38.

To ensure a smoother contour transition, the designed curve was moved downstream, and the arc-circle expanded. The final TIC contour is displayed in Figure 39.

These changes led to some variations in the exit height, nozzle length, and exit angle. A comparable TOP contour was then traced for the design parameters of Table 2.

The performance of both the TIC and TOP nozzles are presented in Table 3. It can be concluded that the TIC contour is very efficient but slightly less than the corresponding TOP contour. This difference can be explained by a smaller average $\theta$ and Mach number for the TOP contour. An oblique shock was also encountered inside the TIC nozzle, as expected, as shown in Figure 40 according to Hoffman [41].

One disadvantage of the TIC design is that the developed algorithm works better for higher Mach numbers and length fraction f. Sometimes, the truncated nozzle is already shorter than the equivalent parabola, but allows for higher divergence at the exit.

As explained before, TIC nozzles can induce critical side-loads, especially during start-up in overexpanding operating conditions.

16.4. Ideal Contour for Wind Tunnel Applications

Wind tunnels are crucial elements in the aerospace industry and research by offering a relatively low-cost mean of simulating real events experienced by vehicles. With the advent of the space exploration era and advancements of aerial aircraft, supersonic and even high-enthalpy hypersonic wind tunnels have become major assets, posing engineering challenges, notably in designing a divergent nozzle capable of providing high quality parallel uniform flow.

Figure 41 illustrates the overall architecture of a supersonic wind tunnel. The primary focus of this Subsection is the design of the divergent D’Laval nozzle, which, like other bell nozzles, typically consists of an expansion section and a straightening section responsible for counteracting deviations caused by the former. In wind tunnel applications, it is crucial that none of the sections have discontinuities, as this would lead to divergences or non-uniformity at the test section.

In the 1930s and 1940s, attempts to trace streamlines through geometry methods were introduced, but proved to be inefficient and often introduced inaccuracies. A faster and more practical analytical method was later presented by Foelsch [27]. Sivells created a sophisticated method formalized in the algorithm CONTUR, already incorporating a boundary-layer correction, also based with the MOC, .

For hypersonic tunnels, where boundary-layer correction and varying $\gamma$ are insufficient to model high-enthalpy flows, the divergent is often designed using Sivells’ method with further optimization through CFD, preferentially involving chemical dissociation modulation.

Potter [42] reported the design of two wind tunnels for testing $M=9$ and $M=10$, utilizing Sivells’ method with some modifications based on the method proposed by Cresci and employing a boundary-layer correction algorithm. However, contrary to expectations, the recorded Mach numbers at the test section were $M=10.15$ and $M=9.3$, respectively.

Benton [43] discovered, through experimental testing, that for high Mach numbers, boundary-layer correction algorithms fail to simulate boundary-layer thickness that has lengths close to the transversal section diameter. For $M=13.5$ and $M=17$, Navier-Stokes equations seem to fail to capture shocks induced by inflection points on the contour, proposing the usage of Parabolized Navier-Stokes Equations.

NASA’s HYPULSE shock tunnel, used to test ramjets, was designed with a code that combines the MOC with an Euler Equations solver, disregarding viscosity, allowing for a considerable test section diameter with uniform and parallel flow [44].

Contemporary wind tunnel nozzle designs often involve the Sivells method, which exhibits good accuracy, especially for high Mach numbers. These designs are complemented by CFD modulation, utilizing sophisticated optimization algorithms such as surrogate-assisted population evolution [21], among others.

16.4.1. Foelsch

In a time when only graphical methods were available, Foelsch [27] proposed a new analytical method to draw the contour of a divergent capable of delivering parallel, uniform flow at the exit section. This method takes its basis in the MOC, although it does not need to be explicitly implemented, but only requires some Equations derived from it. Figure 43 illustrates the philosophy behind this analytical method, whose major goal is to convert the radial flow of an axisymmetric nozzle into parallel flow.

Region TREA is responsible for the expansion of the flow, where it still is radial. As for region AEB, the flow must turn from axial to parallel, as is often called the transition zone. The contour can be traced by projecting left characteristics from points on the right characteristic AE with mass conservation considerations. Since in this region, expansion is negligible, characteristics PK can be considered linear, treated as characteristics in 2D flow.

The drawing of the transition curve, between A and B, is then a simple exercise of fixing several points P, with associated values of ${\theta }_P$, $r_p$, and $M_A\le M_P\le M_E$. The characteristic EB is a straight line with constant Mach number (equal to the desired one for the test section) and parallel flow. For each discrete value of $M_P$, the corresponding point K on the contour is given by Equations 176 and 177.

$$x={\displaystyle \frac{D}{4sin(\omega /2)}}{\displaystyle \frac{\sqrt{ϵ\left(M_P\right)}}{\sqrt{ϵ\left(M_E\right)}}}{\displaystyle \frac{1+\left[cos{\theta }_P\sqrt{M_P-1}-sin{\theta }_P\right]F\left({\theta }_P\right)}{sin{\theta }_B\sqrt{M_P-1}+cos{\theta }_P}}$$

(176)

$$y={\displaystyle \frac{D}{4sin(\omega /2)}}{\displaystyle \frac{\sqrt{ϵ\left(M_P\right)}}{\sqrt{ϵ\left(M_E\right)}}}F\left({\theta }_P\right)$$

(177)

$$F\left({\theta }_P\right)=\sqrt{{sin}^2{\theta }_P+2(cos{\theta }_P-cos\omega )(\sqrt{M_P{}^2-1}sin{\theta }_P+cos{\theta }_P)}$$

(178)

The initial expansion zone can be approximated as an arc-circle, as seen in Figure 43, with radius R, followed by a straight line, with length $\lambda$, until meeting the desired area ratio for point A. To ensure a more user-friendly algorithm, it is set $R=\lambda$ obtained with Equation 179, and the conical angle $\omega$ will be the slope of the straight line. The arc-circle is traced with Equation 180 varying $0\le x\le Rsin\left(\omega \right)$.

$$R=\lambda ={\displaystyle \frac{D}{4\sqrt{ϵ\left(M_E\right)}sin(\omega /2)}}{\displaystyle \frac{\sqrt{ϵ\left(M_A\right)}cos(\omega /2)-1}{cos(\omega /2)+sin(\omega /2)}}$$

(179)

$$y=R_t+R-R\sqrt{1-{\left({\displaystyle \frac{x}{R}}\right)}^2}$$

(180)

Lastly, the two curves must be forced to meet at point A, being brought to the same referential by subtracting from the abscissa of the transition curve a quantity $x_0$ defined in Equation 181.

$$x_0={\displaystyle \frac{D}{2\sqrt{ϵ\left(M_E\right)}tan\omega }}-Rtan(\omega /2)$$

(181)

The typical algorithm for tracing a Foelsch wind tunnel axisymmetric nozzle involves the following steps:

• Have as design parameters the adiabatic index $\gamma$, the test section Mach number $M_E$ and the exit diameter $D\ge 2\sqrt{ϵ\left(M_E\right)}$;
• Calculate ${\nu }_E$ and $\sqrt{{ϵ}_E}$
• Assume $\omega ={\displaystyle \frac{1}{2}}{\nu }_E={\theta }_A={\nu }_A$;
• Calculate $M_A$ and $\sqrt{ϵ\left(M_A\right)}$;
• Determinate the transition curve varying $M_A\le M_P\le M_E$;
• Calculate R, $\lambda$ and $x_0$, subtracting the latter to the transition curve abscissa;
• Trace the arc-circle.

An example of a Foelsch nozzle can be seen in Figure 44. When compared with a minimum length asymmetric nozzle, it is found that the Foelsch method provides a slightly longer one. The length increases with the exit diameter, as expected, but the connection between the two curves gets odd, as seen in Figure 45.

A smoother initial curve can be obtained by changing the position of the beginning of the transition curve, as per Equation 182, and from there, tracing a new initial curve, as in Equation 183. The results of this corrected algorithm are seen for the same nozzle as in Figure 45 in Figure 46.

$$x_A={\displaystyle \frac{3(y_A-R_t)}{2tan\omega }}$$

(182)

$$y=R_t+x^2{\displaystyle \frac{tan\omega }{x_A}}\left(1-{\displaystyle \frac{x}{3x_A}}\right)$$

(183)

16.4.2. Sivells

Foelsch, alongside other authors, proposed methods that focus on building an analytic curve, which was reported to cause inaccuracies responsible for generating shocks that jeopardize parallel uniform flow in the test section. Even if the contour is continuous without major discontinuities, it causes disruptions in the Mach number along the centerline of the nozzle, and induces discontinuities in radial velocity gradients.

Some methods were proposed to reduce the radial velocity gradients responsible for some properties discontinuities, but it was only with the increase in computational power that Sivells [46] proposed a new method. This method involves dividing the axis into three sections, where the Mach number follows a particular polynomial and serves as a boundary condition for the MOC. The contour is then traced with mass conservation equations derived from the MOC.

As seen in Figure 47, a transonic solution is established at the circular throat, in this case, the one reported by Hall [7]. This establishes the sonic line deviation $ϵ$ and leads to the first characteristic, which is the axis reflection of the sonic line, at point I.

From point I to point E, it is assumed that the axial Mach number respects a fourth-degree polynomial. From E to B, the flow is radial, and the axial Mach number has a linear distribution, leading to GA being a straight line with inclination $\omega$ with the axis. The characteristic lines GE and AB respect the assumption of Equations 184 and 185. $\omega$ values should not exceed 15º and are typically around 12º.

$${\nu }_G={\nu }_E+2\omega$$

(184)

$${\nu }_B={\nu }_A+2\omega$$

(185)

Lastly, from B to C, the axial Mach number respects a fifth-degree polynomial, with CD being a straight line with $\nu =\nu \left(M_{exit}\right)$ and $\theta =0^{\circ }$. It’s important to note that the polynomial’s coefficients are chosen so that they coincide at points E and B, and the function has a continuous second derivative that is null at point C.

Instead of an in-house code, Sivells’ FORTRAN code, named CONTUR, as described in Reference [47], was utilized. A Python interface, as a library [48], was used to connect the in-house algorithms for further use of this code.

There’s an inherent limitation in the inputs; for example, only contours for $\gamma =7/5$ can be traced. However, this is not a significant issue as it is a reasonable value for low-density air, typically used in wind tunnels. The contour following Sivells’ method obtained through this code is obviously more reliable and trustworthy.

Figure 48 shows a Sivells contour produced in CONTUR. At first glance, it appears more lengthy than the equivalent Foelsch nozzle but with a smoother contour. Figure 49 compares the same Sivells contour but using the boundary-layer correction supported by CONTUR.

17. Method of Characteristics for Advanced Nozzles Design

Single-stage-to-orbit vehicles (SSTO) are an exemplary concept of space vehicles with simple and efficient system design. But conventional nozzles are unable to provide the highly required efficiency to propel these state-of-the-art vehicles.

Due to their fixed geometry that boundaries the flow, conventional nozzles can’t compensate the varying pressure during the rocket’s ascent. Many solutions have been proposed throughout the years, with plug nozzles, mainly aerospikes, dual bell nozzles, E-D nozzles and Multiple Grid Nozzles (MGN) being the most developed ideas [49].

18. Aerospike

An aerospike nozzle is a type of plug nozzle, where the center-body, or plug, assumes the shape of a "spike" . The term "aerospike" actually refers to the existence of a truncation, where the solid spike is replaced by one that is aerodynamically fashioned [50].

The aerospike nozzle does not have any moving parts and obtains its altitude adaptive properties by lacking a rigid outer boundary, unlike conventional nozzles [51]. This allows it to adjust to the ambient pressure and meet the necessary pressure requirements without any constraints.

The nozzle’s mass is reduced by having a truncated spike, which also helps to avoid some extreme heat fluxes. However, this comes at the cost of the loss of some thrust delivery [49].

During Lockheed Martin X-33 space shuttle development, an aerospike engine (Rocketdyne XRS-2200) was tested onboard a SR-71, with promising results. Eventually this aircraft program was abandoned in 2001 [50]. Since then, some cold and hot tests have been performed, basing most sources of comparison for numerical experiments.

As for the altitude adaptation capabilities of the aerospike, it must be noted that the exhaust gases are not expanded against a bounding wall, as in a conventional nozzle, so besides Equation 1, another term must be added so the force that the gases exert on the spike can be considered, being given by Equation 186. For truncated spikes, the base area also contributes for thrust production, as seen in Equation 187 [52].

$$F_{spike}=\int A_{spike}(p_{spike}-p_{atm})dA$$

(186)

$$F_{base}=(p_{base}-p_{atm})A_{base}$$

(187)

At low altitudes, as shown in Figure 50a), the high ambient pressure causes the exhaust flow to remain close to the nozzle’s wall and be directed to the axial position through a series of expansion-compression fans. The flow is pressed against the wall to the extent that an open wake at ambient pressure is formed. Without truncation, shocks can occur, resulting in flow separation and reattachment, which will lead to undesired side loads.

At higher altitudes, as seen in Figure 50c), the free boundary promotes the appearance of recompression shocks that enable the exhaust flow to remain primarily axial, even in a wider plume. The flow does not expand past the nozzle like in other underexpanded nozzles, thus contributing to thrust production.

At the design altitude, seen of Figure 50b), the wake is already closed, but the major difference from the underexpanded scenario is the plume is less wider, with a more axial flow, promoting a slightly higher thrust.

A secondary gas flow can be introduced into the truncated area to minimize losses arising from the transition from open to close wake [51].

The truncation of the spike also enhances the rocket’s overall performance by reducing the weight of the nozzle. Nevertheless, it should be noted that the transition between open and closed wake occurs at lower altitudes as the truncation magnitude increases. Also, as the truncation is increased, and the base area augments, the advert pressure gradient at the open wake base has its effects amplified according to Equation 187.

At pressure ratios lower than the design point, the recirculation promotes a wall pressure lower than the ambient pressure, resulting in reduced thrust production. However, above the design point, the recirculation zone has a higher pressure than the environment, which positively contributes to thrust production. Therefore, the magnitude of the truncation must be carefully studied to optimize overall performance [49].

Both toroidal and linear aerospike nozzles provide relatively efficient pressure adaptation, with almost ideal behavior below the design altitude. Still, their behavior above the design point is still non-optimal and even performs worse in terms of thrust due to clustering and truncation, compared with bell nozzle designs for higher area ratios [49].

When a wider analysis is performed, at vacuum operation, aerospikes provide a much lighter, shorter and easier to integrate nozzle compared with long bell nozzles. The thrust loss is compensated by providing an improved specific impulse [53].

H. Greer [54] developed a method for rapidly designing an aerospike nozzle contour using the centered expansion fan and geometry relations that relate streamlines with solid surfaces. Basically, this method involves drawing a specific streamline that will serve as the contour of the spike.

The first assumption to implement this simplified method (and often also assumed in traditional MOC) is that at the throat, the flow is sonic in a straight line and not a curved line.

As seen in Figure 51, it must also be assumed that the flow expands solely due to the centered expansion fan on lip A. Mach lines, which are characteristic lines, will expand and straighten the flow so that the desired exit Mach number is achieved with full parallelism to the plug centerline.

The last assumption, more related to the design, is that if the flow is desired to be parallel to the centerline at the exit, the flow at the throat must be incident with an angle equal to the Prandtl-Meyer angle of the chosen exit Mach number.

G. Angelino [55] took Greer’s work and developed an even simpler method where the contour points are calculated discreetly, with an improved accuracy expected. To do so, for an interval of Mach number from 1 to the desired exit Mach number, the length of the respective Mach line from lip A until it meets the streamline that passes through lip B is calculated, as well as the inclination angle in relation to the horizontal.

Angelino’s method has two steps for each Mach number to determine the position of the corresponding contouring point: calculating the distance from lip A to the contour streamline and the respective inclination. The length of such a line for a given Mach number can be obtained in a dimensionless value by the throat length, according to Equation 188, which relates the area of the throat with the area crossed by the velocity timeline that has that given Mach number. It can be further simplified into Equation 189, with $ϵ\left(M\right)$ being given by Equation 30 and $\lambda$ such dimensionless length. Lastly, the inclination is given in Equation 190, which is easily obtained with geometric relations.

$${\displaystyle \frac{l_{line}}{l_t}}={\displaystyle \frac{A}{sin\mu }}{\displaystyle \frac{1}{A_t}}$$

(188)

$$\lambda \left(M\right)=M\times ϵ\left(M\right)$$

(189)

$$\alpha \left(M\right)=\nu \left(M_e\right)-\nu \left(M\right)+\mu \left(M\right)$$

(190)

The simplified method proposed by Angelino provides a nozzle contour that is pretty similar to the one obtained from the full MOC, as demonstrated in Figure 52. The produced contour only deviates slightly by the end of the nozzle, which can be truncated. Additionally, it is observed that the ratio between the pressure experienced by the flow and the stagnation pressure is quite coincident, with the bigger deviations happening at the end of the spike, which in many designs can be truncated.

Nonetheless, Angelino simplified method calculates discreet points that need to be united. More complex methods, such as considering the points immediately after the throat are part of a circumference, don’t produce much different contours then just interpolating all the points. In Reference [56] not only this is proven as also it is concluded Angelino’s simplified method agrees with the results from more complex commercial numerical analyses tools and hot fire test experiments.

More examples of the usage of this simplified method are found in studies by Kumar [57], Dakka [58] and Padania [59].

By inserting as inputs $M_{exit}=4.16$, $\gamma =1.4$ and 100 lines, the designed contour depicted in Figure 53 is obtained. The corresponding data and dimensions of the aerospike are provided in Table 4. It is important to note that the dimensions of the aerospike already take into consideration a throat gap length of 4.005 mm.

The geometric dimensions of the aerospike coincide with those of Reference [60] ($H=2$ in and $L=8$ in), and the contour of the wall, as shown in Figure 54, also exhibits close resemblance. This similarities confirm the accuracy of the Python code used to generate the contour, as it produces the same results as the consulted bibliography.

Numerical studies were further developed with his contour, with all the methodology, results and aerospike altitude adaptative evidences described in Reference [89]. It demonstrates that the aerospike nozzle adapts efficiently to altitude variations, maintaining nearly axial exhaust flow and reducing performance losses compared to conventional nozzles. The study finds that specific impulse increases with altitude, with minimal losses in underexpanded conditions but noticeable reductions in overexpanded scenarios due to turbulent effects. The realizable $k-ϵ$ turbulence model is validated as a reliable tool for simulating aerospike flow characteristics, reinforcing the nozzle’s advantages for space propulsion applications.

In summary, the advantages of the aerospike nozzle include:

• Vector Thrust - Due to clustering of at least four different throats, it is possible to induce pitch, yawn and roll (just for linear) moments. More complex methods, such as flapping, are necessary to induce roll moments in toroidal aerospikes [61]. This eliminates the need for a gimbal or other movable piece, reducing the weight of the engine [62].
• More Average Impulse - At low altitude, the aerospike behaves almost like the ideal nozzle, compensating more substantial losses felt at higher altitudes. This proves aerospike nozzles are suitable and advantageous for low altitude operation and also for SSTO vehicles [63].
• Better Rocket Aerodynamics - The base of the rocket produces drag. In an aerospike, almost all this area is used to produce thrust instead. Since all the base can be used to produce thrust, at vacuum operation, area ratios of up to 1:150 can be obtained [63].
• More Efficient Rocket Structure - The thrust produced by the nozzle is transferred to the rocket as a distributed load and not as a point source load, allowing for a more efficient and lighter rocket’s structure [63].
• CubeSat application - Aerospikes can be design for vacuum operation, being easily scaled and integrated in CubeSat designs. They will be much shorter than a bell nozzle with the same area ratio. Vector thrust, high area ratio, and low weight make a good trade-off for using aerospike engines in satellites with a low form factor [53].
• Better Packing - When multiple staging cannot be avoided, aerospikes are much easier to transport to high altitudes as for high expansion ratios, have much shorter length and volume compared with an equivalent bell nozzle [64].
• Improved Reliability - Clustering allows the several modules to operate independently. In case of a malfunctioning module, it can be shut down as so an 180º opposing one.This way, the mission is not jeopardized, keeping the thrust vectoring with just the loss of some thrust [64].
• Practical Landing - Clustering allows the several modules to operate independently. During landing, a reduction in thrust can be obtained by sequentially shutting down pairs of modules or reducing the flow in some. This allows a more reliable landing method with only the need to develop controlling valvules. There’s a reduction in “legs” weight as they can be shorter due to the shorter nozzle [64].

As for the major disadvantages and possible was to overcome them:

• Extreme Flow Condition - The truncation experiences high temperature, pressure, and heat fluxes, and its small volume makes it difficult to implement regenerative cooling and other solutions. The usage of $90\%$ H2O2 fuel leads to a decrease in the temperature of the exhaust gases to around 1023K, which can be easily withstood by steel alloys [65].
• Considerable weight (full spike) - The solid spike is typically made of heavy metal elements with complex coatings that can withstand the extreme conditions of the expanding flow. Truncated nozzles can be made shorter, by truncation, but there is a trade-off between weight and average impulse [63].
• Negative Base Pressure - Closing the wake before the design point, results in a lower base pressure than the ambient pressure, which reduces thrust, as explained previously. secondary flow, prevenient from the gas-generator-cycle that powers the tanks pumps can be added at the base, increasing the local pressure and, therefore, the thrust produced [66].
• Clustering Flow Interaction - When the flow from each module expands and interacts, it creates additional shocks, which decreases the produced thrust. No solution was found but the performance losses are small and does not mean bell nozzles are more efficient [67].
• Manufacturing challenges - The small area of the throat, coupled with the extreme heat flow, makes it difficult to assure continuous optimal performance [68]. Also, in most designs, the spike needs to be fixed inside the combustion chamber. Additive manufacturing techniques can produce alloys and other advanced materials that are capable of withstanding the harsh environment around the nozzle, while maintaining the desired dimensional accuracy [69].

19. Dual-Bell

First proposed in 1949 by Cowles and Foster [70], dual bell nozzles are a very interesting concept that delivers altitude adaptation behavior by varying the effective expansion ratio of the flow, discretely. Depicted in Figure 55, this concept consists of two bell nozzles, with the first designed for low altitudes, followed by an attached extension that provides higher expansion for vacuum operation.

In a mission proposed by Immich [71], comparable to that of an Ariane 5 equipped with a dual bell nozzle, a $72\%$ payload increment, from an original 2000 kg, can be obtained. Ferrero [72] finds a 1500kg payload increment, already taking into account modern missions and contemporary understanding of the dual bell nozzle’s flow behavior. Taylor [73] showed that the dual bell is preferable to the E-D nozzle for overexpanding scenarios, with competitive overall total impulse.

The dual bell nozzle offers two modes of operation [49]:

• Mode 1 (low altitude operation) - Due to high wall pressure at the extension, aligned with the overextended flow from the inner nozzle, this mode allows for symmetrical flow separation at the inflection between the two bell contours. This way, the effective expansive area ratio of the nozzle is smaller, reducing overexpansion effects and related losses.
• Mode 2 (vacuum operation) - Lower wall pressures at the extension and flow’s underexpansion make it adhere to the extension, fully utilizing all the nozzle’s effective expansion area, thus reducing specific impulse losses.

In theory, the performance of a dual bell nozzle should be superior to the usage of a single bell nozzle with the equivalent area ratio of each segment. Analyzing the specific impulse curve of the dual bell nozzle throughout an ascending flight envelope, in Figure 56, reveals some additional losses.

In mode 2, the specific impulse is slightly inferior to a Rao nozzle with the same exit divergence and area ratio, due to the imperfection in contour caused by the inflection point that induces a weak oblique shock, as seen in Figure 57 b).

Losses of similar magnitude occur in mode 1, where the align effect of the aerodynamic performance of the rocket, which decreases the pressure at its base, and the overexpansion of the flow from the inner nozzle decrease the wall pressure of the extension, as indicated by Figure 57 a). The effect of this aspiration drag decreases during the ascent of the rocket but can be up to $6\%$ at ground level.

The most significant loss of specific impulse results from an early transition between modes, provoked by the aspiration drag that reduces the wall pressure of the extension. The transition itself can lead to other engineering problems, such as side-loads [49].

Recent studies have focused on both delaying and accelerating the transition between modes to minimize the impact of side-loads [71]. Constructing structures capable of resisting propagated side loads is not recommended, as it would consume most of the payload gains and introduce flight stability challenges.

At the inflection point, the flow is bent outwards, promoting its expansion, which can further decrease the wall pressure. Therefore, the extension section must be designed with consideration for the desired wall pressure gradient it needs to promote [74]:

• Favorable $\Delta P<0$ - Similar to whats observed in traditional bell nozzles, decreasing pressures will lead to uncontrolled attachment and reattachment.
• Null $\Delta P=0$ - In theory, it promotes an instantaneous transition as soon as the wall pressure, as a whole, drops below the critical value.
• Favorable $\Delta P>0$ - Also promotes an instantaneous transition, but in traditional nozzles, it generates a tremendous side load, which is hardly felt if the transition occurs fast enough.

Taylor [73] proposed achieving a later transition by using an inner nozzle projected for higher Mach numbers, acknowledging a compromise with additional losses at low altitudes. Experimental studies demonstrated that a 0.4-second transition could be achieved without registering side loads, although exit velocity profiles showed evidence of asymmetrical displacement.

Adjusting the flow’s pressure can be a method to promote a later transition almost instantaneously. Frey [74] demonstrated that chamber throatily can be employed by constantly decreasing stagnation pressure throughout the rocket’s ascent and then increasing it back to the original value instantaneously when the transition is desired. Ferrero [72] suggested fluidic control of the transition by adding high-pressure flow near the inflection point, creating a barrier. This mechanism would be activated near the premature transition altitude and shut down once full contour utilization is desired. Numerical studies confirmed that this method provides a perfectly timed transition with minimal registered side loads.

Wang [75] proposed a combination of a dual bell nozzle with a central pintle typical of an E-D nozzle. In overexpanding scenarios, this combination performed better than the E-D nozzle but worse than a dual bell. The surface of the pintle, having low pressure, accentuated aspiration drag effects, although it was compensated with a later transition. Experimental findings revealed that CFD models underestimate specific impulse, as the simulated wall pressure was lower than experimental values. Additionally, in underexpansion, near the exit lip, some zones with favorable pressure gradients were observed, which could jeopardize a perfect transition.

Kbab [76] and [77], and Wu [78] presented additional numerical studies that confirmed the previous concepts, findings, and assumptions. Moreover, it was noted that the MOC tended to over-predict the expansion at the inflection, and boundary layer correction to the contour provided more accurate results [76].

In terms of design, literature suggests that both TIC and TOP nozzles can be used for the contours of each segment of the dual bell. If a specific pressure gradient is desired, the extension can be designed with the MOC. For example, Taylor [73] and Kbab [76] use the MOC to trace the isobaric line emanating from the inflection point, which will serve as the extension contour.

An in-house code was developed for designing a dual nozzle in the same fashion as Otsu [79], but with an arc-circle downstream of the geometrical throat, with ${\theta }_{N1}={30}^{\circ }$, as it is a typical value for traditional bell nozzles.. Depicted in Figure 58, it consists of two TOP nozzles.

Both nozzles are referenced to a conical nozzle with a half angle of $\alpha ={15}^{\circ }$, a length corresponding to $f=0.8$, and accepting an exit divergence angle of ${\theta }_{E1}={\theta }_{E2}={15}^{\circ }$.

The inner nozzle is designed for an expansion area ratio of 37.5, and the outer nozzle for 100. The extension initial angle is given by Equation 191

$${\theta }_{N2}={\theta }_{E1}+\beta ({\nu }_2-{\nu }_1)$$

(191)

Otsu [79] investigated the acceptable values of $\beta$, which controls the extension’s wall pressure gradient by modulating the Prandtl-Meyer angles difference at the exit of the extension and inner nozzles. It was found that values of $\beta \le 1$ lead to long transitions between modes, starting at higher ambient pressures. For $\beta =1.2$ and $\beta =1.5$, not only was the transition fast, taking only 3 ms, but it also occurred at a lower ambient pressure, closer to the ideal value.

The depicted nozzle has $\beta =1.2$ and ${\theta }_{N2}=26.{05}^{\circ }$, representing a good trade-off between a swift transition and manufacturing feasibility.

In summary, the advantages of the dual bell nozzle include:

• Doesn’t have movable parts, allowing for traditional cooling systems, increasing reliability [49].
• Payload gains [71,72].
• Net impulse gains in comparison with reference traditional Rao nozzles [49].
• It can be implemented with any traditional combustion chamber [73].
• Higher thrust to engine weight ratio [80]

As for the major disadvantages:

• Low rocket’s base pressure, creating aspiration drag [49].
• Premature transition can degrade the specific impulse of the total envelope [71].
• Dangerous side loads can be generated during transition [74].
• At vacuum operation it behaves as any high area ration traditional nozzle but with additional losses [73].

20. Expansion-Deflection

Rao [81] proposed the concept of the E-D nozzle, as depicted in Figure 59. It consists of an annular axisymmetrical throat that is tilted towards a contoured wall. The gases are expanded around the central plug or pintle, and the contour ensures they maintain a mostly axial direction at the exit plane, where flow conditions determine the thrust produced.

The name E-D stands for the expansion promoted by the central pintle and the deflection imparted by the contour. It operates in an altitude-adaptive manner similar to that of plug nozzles, with the difference being that the expansion is controlled from within the nozzle, as the throat delivers flow outwardly, meaning there’s an objective area ratio ceiling.

Figure 60 illustrates the flow inside an E-D nozzle, where a system of compression and decompression waves directs and adapts the exit pressure of the flow. An open wake at the center of the exit flow creates aspiration drag, reducing the pressure in the wake to which the flow pressure will equal, leading to overexpansion.

As the outside pressure drops, there’s an increase in the area ratio, closing the wake, or, more precisely, the pressure of the expanding flow is now higher than that on the wall of the base of the pintle. However, the losses incurred from overexpansion negatively impact the altitude adaptation of E-D nozzles, making them only more efficient when compared to high area ratio bell nozzles with a desired short length [49].

Since, once the wake is closed, the gases at the base of the pintle stay trapped and create some recirculation flow, eventually contributes positively to thrust. Taylor [82] argues that almost vacuum performance can be achieved with a short E-D nozzle with the contour designed using inviscid methods for optimum thrust using the MOC. It is also suggested that the open wake mode is characterized by high levels of turbulence and other isentropic losses.

The nearly optimal behavior of E-D nozzles at vacuum operation, characterized by a relatively small area ratio, is later demonstrated through cold testing by the University of Bristol [73].

The high viscosity effects felt during open wake are confirmed by Schomberg [83], who concludes that contemporary CFD methods are accurate to simulate E-D nozzle flow.

The central pintle is a critical aspect of E-D nozzles, significantly influencing their performance by defining the transition point and its process. Paul numerically experimented with some pintle geometries by contouring their expanding surface, such as an arc-circle with several radii. For open wake operation, there’s little influence as there’s a separation point right after the throat. A higher radius will lead to a sooner and softer wake closing, which is advantageous as aspiration drag will start to decrease. Nonetheless, it is only at even lower ambient pressures that the pintle base will start to contribute positively [84].

Taking into account literature suggestions, no E-D nozzle was projected, especially for full flight envelope operation, as there are many design considerations with no inviscid method since both the pintle and contour need to be carefully studied.

In summary, the advantages of the E-D nozzle include:

• Doesn’t have movable parts increasing reliability [49].
• It can produce the same thrust of a traditional bell nozzle but with less than half the length [81].
• Structural advantages if the throat is near the axisymmetric line [81].
• Excellent in vacuum operation [82].
• Compared tot he aerospike nozzle, it is more easily installed in traditional combustions chambers [82] .

As for the major disadvantages:

• High heat flow around the pintle [7].
• Low altitude adaptation capabilities due to aspiration drag from the open wake [82].

21. Multiple Grid Nozzle

With a transversal sketch depicted in Figure 61, where each triangle represents a "nozzlettes", MNG have been studied in recent years, not necessarily for their altitude adaptive behavior, but also for their potential to improve overall performance.

The way a MNG increases the expansion ratio is by adding more "nozzlettes" rather than increasing the length, which adds little weight and has been shown to improve performance when designed with multidisciplinary methods.

Anyhow, the small throat area of each "nozzlettes” leads to high ablative erosion, specially to the throat, which will drastically decrease performance over time [85].

22. Other Advanced Designs

Some nozzle designs have been proposed throughout the years but with little contemporary research.

22.1. Nozzles with Fixed Inserts

This type of bell achieves a similar behaviour and performance curve as the dual-bell nozzle by inserting a trip ring inside a conventional bell nozzle, responsible for causing the flow separation at low altitudes.

According to some literature [86], this ring can be as small as 10% of the local boundary-layer thickness in order to reduce the losses at vacuum performance.

The transition point between the sea level and vacuum operation depends on the wall pressure near the trip ring and the disturbance provoked by it, creating significant uncertainties. Although this is a simple and with low technology risk concept, the referred uncertainties led to the rejection of this concept [49].

22.2. Nozzles with Temporary Inserts

The presence of an insert within the nozzle negatively impacts vacuum performance. If the insert is removed before high altitudes are reached, the high performance of the baseline nozzle can be matched.

A modified RD-0120 engine with a secondary nozzle insert was previously hot-fire tested, resulting in a 12% performance gain at sea level, with the only losses compared to the ideal nozzle being caused by aspiration drag.

The challenge with this design is the ejection of the insert, which can lead to undesired side loads, shocks, and collisions near the exit plane. Although ablative inserts have been considered, the stability of the surface and the regression rate are still difficult to model with current technology [87].

22.3. Nozzles with Active Secondary Gas Injection

This concept proposes that a second gas flow can be injected inside the nozzle, normal to the wall, in order to force flow separation during sea level operation, reducing the nozzle’s area ratio.

The injection is originated by a high-pressure reservoir, but the amount of fluid needed to induce flow separation is considerable, without expected net impulse gains [49].

22.4. Nozzles with Passive Secondary Gas Injection

Similar to the previous concept, but the injection is achieved by creating “slots” to let the ambient air enter the nozzle. The "slots" are closed at high altitudes.

Unfortunately, this passive mechanism only works at lower altitudes, where the ambient pressure is high enough to enter the nozzle. This operational envelop is further shorted due to the stream around the rocket, that reduces the pressure at its end, near the engine’s nozzle [49].

22.5. Extendable Nozzles

This concept, seen in Figure 62 involves a truncated low area ratio nozzle that is used in sea-level conditions. It has an outer extendable component that attaches at higher altitudes, thereby increasing the area ratio. The nozzle can be design using the Method of Characteristics, for example, showing an improved performance.

However, the in-flight placement of the extendable component requires the use of mechanical devices, which increases the nozzle’s weight and reduces its reliability, posing engineering challenges for the cooling system. Some upper stages of rockets already use a similar concept to reduce packaging volume.

22.6. Nozzles with Throat Area Varied by a Mechanical Pintle

The mechanical pintle, used to control thrust in some solid propellant rockets equipped with bell nozzles, is a proposed method to achieve altitude adaptiveness, for this and other types of propellant. The pintle works by adjusting the throat area, which in turn would alter the nozzle’s area ratio.

However, the implementation of such a mechanism requires a complex and heavy actuator and cooling system, leading to increased weight and decreased reliability [49].

22.7. Dual-Throat Nozzle

The concept of the dual-throat nozzle involves two conventional combustion chambers, one concentric to the other, with both converging to the same fixed exit area, as depicted in Figure 63. The original purpose of the dual-throat nozzle was to regulate thrust, but it can also be altitude-adaptive.

At low altitudes, both combustion chambers are active, resulting in a larger throat area and a lower area ratio, thereby preventing overexpansion. At a certain altitude, the outer chamber is shut off, reducing the throat area, increasing the area ratio, and allowing the flow to fully expand.

Both behaviors can be seen in Figure 64. Still, there are some losses in vacuum operation due to the discontinuity caused by the outer chamber and the resulting recirculation [49].

22.8. Dual-Expander Nozzles

Similar to the dual-throat nozzle concept, the dual-expander nozzle has the same geometry and altitude compensation behavior, with the difference that the inner chamber has a short divergent section and is the one that is shut down during vacuum operation.

Both behaviors can be seen in Figure 65. It is important to highlight the rapid expansion and subsequent recompression in Figure 65b), as it creates a recirculation zone. This issue can be addressed by bleeding air into the inner chamber, though there will be a small reduction in impulse due to the increased area ratio.

The major engineering challenges associated with this design include the potential structural disintegration of the inner chamber during vacuum operation due to high heat fluxes and pressure oscillations [49].

23. Conclusion

This research has advanced the understanding of traditional and modern nozzle designs, contributing valuable insights into their operational efficiencies under varied flow conditions. The investigation of conical and bell nozzles in Chapter 3 demonstrates the effectiveness of the Method of Characteristics in optimizing nozzle contours for different performance metrics, including thrust, specific impulse, and flow stability. Results highlight that parabolic and truncated configurations can achieve near-optimal performance while significantly reducing material and fabrication requirements.

The advanced nozzle configurations presented in Chapter 4 exhibit superior adaptability to varied ambient pressure conditions, with designs such as the aerospike and dual-bell nozzles showing notable advantages in both thrust optimization and flow stability. These designs, though complex in their geometrical requirements, hold promise for future propulsion systems requiring versatility and efficiency. For example, the aerospike nozzle’s altitude compensation and the dual-bell nozzle’s ability to switch operational modes demonstrate innovative solutions for adapting to expanding performance needs in aerospace engineering.

Besides validation and comparison with commercial CFD tools, future research should explore the integration of these advanced designs with new materials and manufacturing technologies, such as additive manufacturing, to further reduce weight and improve structural integrity. Additionally, extending these investigations into real-world experimental validations could bridge the gap between theory and practice, accelerating their adoption in commercial and military aerospace applications.