Contemporary Importance of the Method of Characteristics for Supersonic Nozzle Design

Rocket technology has a long history, initially linked to warfare. Early examples include 1^st-century Chinese bam-boo tubes filled with gunpowder, evolving into “arrows of flying fire”, that in 1232 were used by the Chinese army against Mongol invaders. Significant advancements occurred in the 20^th century, driven by World War II and the Cold War. American engineer Robert H Goddard, launched 34 liquid propellant rockets between 1926 and 1941, reaching altitudes of 2.6km and speeds of 885km/h. He later developed single and multi-stage rocket engines. German engineer Wernher von Braun developed the V-2 rocket for the Axis forces during World War II. After the war, he joined NASA and played a key role in developing the Saturn V launcher that would propel Apolo 11 to the moon in 1969 [1]. Ever since, the optimization of rocket systems, by increasing efficiency and readability, have been pursue and a major concern of aerospace industry. Thrust in aerial vehicles is generated by mass ejection, explained by Newton’s third law and quantified by Equation 1. Rockets have supersonic nozzles installed composed by a convergent section (𝑀<1), a athroat (𝑀=1) and a divergent section (M>1), all responsible for converting the thermal energy of the flow debited by the combustion chamber into kinetic energy, reducing temperature and pressure [2]. Some important nozzle performance criteria are the specific impulse that measures the thrust produced per weight flow rate of propellant burned (Equation 2), the total impulse signifying the thrust produced over the entire engine’s operation duration (Equation 3) and the thrust coefficient which is a dimensionless thrust value obtained as its ratio with the chamber pressure and the throat area (Equation 4) [3].

Nozzle design importance

Weight savings and high reliability are crucial in any rocket design, and the respective nozzles are no exception. The goal of their design is to expand the exhaust flow so it produces the most thrust with the shortest length possible. The ideal nozzle represents the best performance possibly achieved by certain chamber conditions and ambient pressure, with assumptions like chemical equilibrium, perfect gas law, adiabatic and isentropic flow, negligible boundary layer effects, steady flow parallel and uniform at the exit.

The design methods of real nozzles aim to match the ideal nozzle’s performance, but face losses due to real gas phenomena, such as [2]:
A. Non-Uni-dimensional flow-Supersonic flow expands outward and thrust is just effective in the direction to which the vehicle is moving. The contour of the nozzle can straight the flow to a parallel positioning at the exit but increasing the nozzle’s length and inducing isentropic losses.
B. Boundary Layer (BL)-The effective expanding area ratio can be reduced, especially for smaller and high enthalpy nozzles, with a percentage of the exit flow between 2% and 25% being subsonic.
C. Reactive flow-The expanding gases have great changes in their temperature, with molecular dissociation and even incomplete combustion inside the chamber. The isentropic expansion and steady combustion utilized in the ideal nozzle can over predict thrust by 0.5%.
D. Multi-phase flow-Liquid droplets and solid particles can be present in the flow and require acceleration, though they contribute minimally to thermal energy. When these particles are larger than 0.015mm and constitute more than 6% of the flow’s mass, specific impulse losses can increase by 10% to 20%. While less critical, in high area ratio nozzles, extreme and fast expansion can cause some exhaust constituents to precipitate.
E. Nozzle erosion-The aggressive expanding flow can erode the walls of the nozzle, changing its contour and effective area ratio, representing losses of specific impulse of up to 0.7% due to chamber pressure drop.

Besides, nozzles often operate outside the design point, meaning the pressure of the exhaust flow at the exit won’t match ambient one, having a different Mach number and exit angle different from design, increasing the expected losses of up to 15% [4]. Contemporary sophisticated methods are able to precisely predict the refereed phenomena and provide flow field within a nozzle although with expensive computational cost that jeopardize the assessment of multiple design in pre-liminary project phases. Depending upon the chamber’s enthalpy conditions and exit velocities aimed, some real gas effects can be overlooked and easier models implemented.

Transonic solutions

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. It became of utmost importance to develop a sufficiently accurate method to model the velocity vectors of the flow in the throat region, considering the transonic nature of the local flow. The first modeling of the throat region was proposed by Hall [5]. It consists of an expansion in order of 𝑅𝑐^-1, where 𝑅𝑐 is the radius of the contour immediately downstream of the geometric throat. By expanding the series of velocity magnitude to 𝑅𝑐^-3, parabolic, hyperbolic, and circular arc throats are taken into account. Other shapes can implement this method with some adaptations. Hall’s method demonstrated considerable accuracy compared to experimental results for 𝛾=1.4, 𝑅𝑐=5 and two-dimensional flow. The series expansion to 𝑅𝑐^-1 yielded an error of 0.5% and when fully expanded to 𝑅𝑐^-3, the error further decreased to only 0.02%. For 𝑅𝑐<1.5, Hall’s solution starts to deteriorate and for 𝑅𝑐<1, which is common among most rocket nozzles, the method diverges in an oscillatory manner. Some inaccuracies can also be detected in the terms of the third degree of the expansion for 𝑢3 and 𝑣₃. To address these problems, Kliegel [6] presents a similar method, but the series are expansions of (𝑅𝑐+1)^-1. This way, even 𝑅𝑐=0, corresponding to a sharp corner, can be considered. Hall’s solution employs cylindrical coordinates, but Kliegel defends that toroidal coordinates shall 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 𝑅𝑐=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 cylindrical coordinates.

Dutton [7] proposes a (𝑅𝑐+𝜂)-1 expansion, where 𝜂 is arbitrary, respecting the differential equations of motion for cylindrical coordinates when bringing the control volume back from toroidal. In fact, as seen in Figure 1(a), Dutton’s method is the most accurate for 𝑅𝑐=0.625 compared to experimental results. An algorithm with Dutton’s method was implemented. Its sonic line for 𝑅𝑐=0.625 was traced and can be seen in Figure 1(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.

Figure 1:Transonic Line with (a) several methods [7] and (b) In-House algorithm with Dutton’s method.

Main supersonic flow models

From the simplest unidimensional isentropic expansions and compressions to the consideration of a wide plenitude of 3D real gas effects, there are many ways to model supersonic flows. Depending upon the context and desired precision of the simulation, a fitting model must be chosen.

Quasi 1-D flow: Supersonic nozzles expand the gas by providing an area variation to the flow. Similar to a unidimensional flow, a quasi1-D model will also allow for area variations allowing the understanding of the function of each nozzle’s section. The flow’s properties are still considered one-dimensional, meaning they are constant within any cross-sectional section. The nozzle’s contours won’t create shocks as they are assuming to have smooth wall angles, maintaining adiabatic and isentropic behavior. This leads to the relation between area and Mach number stated in Equation 5, suggesting that for 𝑀<1, flow accelerates with decreasing areas, as in the convergent section of a supersonic nozzle. At a 𝑀>1 regime, velocity will increase with increasing areas, as observed in the divergent section. As for sonic velocities, suggests that a minimum area is reached, which in a rocket nozzle represents the throat. When 𝑀=1, 𝑑𝐴/𝐴=0, meaning sonic velocities are reached at a minimum area section, corresponding to the throat of the nozzle.

The area ratio offered by a duct can lead the flow to a supersonic or subsonic Mach number, as peer Equation 6. Flow’s properties at any point in the duct can be modeled considering isentropic expansions [8].

Method of characteristics: Proper nozzle’s contour analysis need to account for two-dimensional effects, so a more sophisticated method than the quasi-1D is needed for proper nozzle’s design and performance assessment. Within a flow, the full-velocity potential in Equation 7, obtained with the combination of the continuity equation and Euler’s equation for a two-dimensional and irrotational flow, takes the form of a Partial Differential Equation (PDE), which is nonlinear, becoming hyperbolic if Equation 8 is satisfied, a condition met in supersonic flow. 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 constant relationship, facilitating the solving of the flow at different points, starting from transonic solutions. To apply the following description of the MOC, the flow must be assumed to be supersonic, inviscid (with neglected BL effects), steady and irrotational.

The direction lines defining the hyperbola are known as characteristics, and, in the context of the velocity potential in supersonic flows, are coincident with Mach lines. As a function of both 𝑥 and 𝑦, the velocity potential’s derivatives are expressed by Equations 9 and 10.

By combining Equations 7,9 and 10, a linear system with variables Φ_𝑥𝑥, Φ_𝑦𝑦 and Φ_𝑥𝑦 is formed. Solving for Φ_𝑥𝑦, using Crammer’s, Equation 11 can be obtained.

At a given point in the flow, Φ𝑥𝑦 implies that along a direction 𝑑𝑦/𝑑𝑥, the velocity undergoes changes corresponding to the respective 𝑑𝑢 and 𝑑𝑣 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 onward a characteristic line, crossing several streamlines. Making the denominator of Equation 11 zero, results in Equation 12. The slope of the characteristic lines that pass through a given point can then be obtained using Equation 13 [9]. Considering Equation 13, 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, an imaginary Mach angle arises, signifying that the characteristic lines take an elliptical form [8]. Examining Figure 2, two characteristics are observed passing through a designated point A: A left characteristic 𝐶₊ with a slope of 𝑡𝑎𝑛(𝜃+𝜇) and a right characteristic 𝐶_- with a slope of 𝑡𝑎𝑛(𝜃-𝜇).

Figure 2:Characteristic lines of point A.

To solve the flow along a characteristic line, the governing PDE describing the flow reduces to ODEs. These equations are obtained by setting the determinant of the numerator in Equation 11 to zero, resulting in Equation 14.

Equation 15 represents the set of ODEs obtained from Equation 14. With algebraic manipulation it culminates to the expressions presented in Equations 16 and 17, for a planar flow.

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

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 20 and 21. In the axisymmetric MOC, the assumption of a constant relation along a given characteristic line no longer holds. Consequently, the relationship between the variables 𝜃 and 𝜈 requires integration, and this is achieved by employing a finite differences method described by Zebbiche [10] and Restrepo [11].

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 and average as depicted in Equation 22. For an axisymmetric case, both the nodes’ properties and their coordinates are solved at once, as there’s no longer a constant through a characteristic line. If the MOC is utilized 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 and most often needs a transonic solution.

Commercial CFD tools: The MOC proves to be a reliable tool for designing supersonic nozzles, mainly due to its simplicity. By considering an inviscid flow, the MOC neglects real gas effects like changing specific heats and the BL. Considering local vibrational equilibrium or applying a displacement to the contour representing the BL are simply implemented methods to improve the MOC. However, they still lack important features specially for hypersonic flow, characterized by high enthalpy effects and where a parallel uniform exit flow is mandatory. Computational Fluid Dynamics (CFD) is a branch of fluid mechanics, including any numerical analysis and 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 real-life prototypes.

Differential conservation equations for inviscid flows: Traditional fluid analysis, and consequently conservation equations, are set for a control volume 𝒱 delimited by surface area 𝑆. So, for any vectorial function 𝐴 and scalar function Φ, Equations 23 and 24, respectively, are true for any control volume.

Equation 25 represents the continuity equation for a given control volume. Taking into account Equation 23 wit 𝐴=𝜌𝑉, Equation 26 is obtained.

For Equation 26 to be zero, one could argue that the first part of the integral could be annulled by the second part. But the analysis must be true for any random control volume, so the only possibility for it to be zero is if the control volume is a single point, precisely what’s desired in CFD algorithms. Equation 27 is obtained, representing the differential equation of continuity.

Equation 28 represents the momentum conservation for a given control volume. If combined with Equation 24 with Φ=𝑝, Equation 29 is obtained. Equation 30 represents its scalar decomposition in the x direction.

Now, using Equation 23 with 𝐴=𝜌𝑉, Equation 31 is netted. Further combination with Equation 30 leads to a single null integral, as depicted in Equation 32.

Solving Equation 32 with the same philosophy of continuity, Equations 33, 34 and 35 represent differential equations of momentum conservation for each Cartesian direction.

Equation 36 is the energy equation for a given control volume. If combined with Equation 23 with 𝐴=𝜌(𝑒+𝑉²/2)𝑉 and 𝐴=𝑒𝑉, Equation 37 is reached.

Solving Equation 37 with the previous philosophy, Equation 38 represents the differential equations of energy conservation [8].

Turbulence models: Typical simulations of the flow field inside a nozzle focus on properly modeling the flow and BL 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 simulations fail, and Reynolds-Averaged Navier-Stokes (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 flow field. Two-equation models, based on Boussinesq’s concept of viscous turbulence 𝜇𝑡, are commonly adopted. Most models include equations for turbulent kinetic energy (k) and either k dissipation (𝜖) or k specific dissipation rate (𝜔).

One widely used turbulence model for nozzle flow is the Shear Stress Transport (SST) k-𝜔, which calculates turbulent viscosity considering the transport effects of the main shear stress. Its greatest advantage lies in gradually transitioning within the BL from a traditional k-𝜔 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-𝜔 model has proven effective for ducted flows from injectors, accurately predicting shock waves and BL detachment, albeit with slight overprediction of shear stresses, potentially leading to higher velocities in secondary flow. Detailed mathematical modeling of SST k-𝜔 can be found in Reference [12] and examples of the utilization in nozzle related numerical studies described by Jia [13], Zhang [14], Lee [15] and Restrepo [11], for both steady and transient flow.

As stated before, the convergent section of supersonic nozzle is responsible to deliver a straight choked sonic to the throat, with no major design concern. Is the divergent section that is most crucial aspect of nozzle design and to which a plenitude of methods has been proposed.

Cone nozzle

The cone nozzle, shown in Figure 3, is the simplest divergent design, being characterized by its area ratio and half angle angle 𝛼. Due to its easy manufacturing, it has historically been one of the most used contours [3]. The design options of the half-angle represent a trade-off between length and divergence angle at the exit. The susceptibility of the cone nozzle for flow separation and overexpansion cannot be overlooked. Besides, isentropic thrust losses increase with the decrease of the divergence angle. When factors such as performance, weight and length are considered, a half-angle of 15º offers the best compromise, with it not surpassing 25º degrees as separation becomes too intense for practical applications [16]. The emergence of CFD methods for nozzle analysis has left a gap in the past decades regarding the application of the MOC to solve cone nozzles, so an in-house MOC algorithm was implemented to gain a deeper understanding of certain phenomena associated. 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. It’s assumed that the flow inside a cone nozzle behaves in accordance with Figure 4, where a centered expansion fan at the lip of the throat is responsible for expanding the flow from 𝜃=0º to 𝜃=𝛼. If the flow at point C has 𝜃=𝛼 and 𝑀=𝑀_{𝑒𝑥𝑖𝑡}, its left characteristic with 𝐾₊=𝛼-𝜈(𝑀_{𝑒𝑥𝑖𝑡}) can be traced until it meets the axis at point B, where 𝐾₊=0-𝜈_𝑏, leading to 𝜈_𝑏=𝜈(𝑀_{𝑒𝑥𝑖𝑡}). From point B, the right characteristic 𝐾_-=0+𝜈(𝑀_{𝑒𝑥𝑖𝑡}) is traced back to the centered expansion fan at point A, where, due to the nature of the referred structure, 𝜈_𝑎=𝜃_𝑎. 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 39.

Figure 3:Cone nozzle’s configuration.

Figure 4:Assumption of the flow inside a cone nozzle.

Unfortunately, the algorithm failed to capture the desired behavior. Firstly, the conditions at point C were observed inside the nozzle, as depicted in Figure 5. 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 𝑀_{𝑒𝑥𝑖𝑡}=2.4 and 𝛾=7/5, the average value between the Mach number at the exit lip and centerline 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 𝛼 is smaller than the one calculated in Equation 39, 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 𝛼. In the case of 𝛼 larger than the one provided by Equation 39, the initial expansion may be excessive, and the characteristic lines from it seem to influence the flow outside the nozzle. At the exit plane a shock may occur in order to raise the pressure of a possibly overexpanded flow.

Figure 5:Characteristic lines inside a 2D ramp designed for 𝑀_{𝑒𝑥𝑖𝑡}=2.4 and 𝛾=7/5.

In reality, cone 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, consistently causing an oblique shock inside the nozzle. This creates an observable pattern of a double Mach diamond, with its thickness peaking near design operation and a small Mach disc forming inside the nozzle as in Figure 6. Shadowgraph images reveal that under expanding conditions cause an increase in expansion at the throat, potentially overtaking 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 [17]. Similar to other traditional nozzles, in overexpansion, the entire contour is not fully utilized, risking the generation of side loads due to uneven separation, especially during start-up [13] and reattachment during the transition to total contour utilization. Historically, to mitigate the presence of the throat shock and achieve an isentropic cone nozzle, the suggestion of incorporating an initial expanding arc has been made. Yet, utilizing the asymmetrical MOC reveals that a weak oblique shock will consistently emanate from the junction between the arc and the cone unless the junction is optimized to render the shock almost imperceptible [18].

Figure 6:CFD Depiction of the flow inside cone nozzle operating at design conditions [17].

Bell nozzles for propulsive applications

Bell nozzles are shorter, with lengths as fractions of cone nozzles, providing smaller contour angles at the exit. These design changes allow a more parallel flow at the exit, improving the overall performance of the nozzle. The gradual reverse of the divergence angle in the shorter nozzle is possible with a 20º-50º divergent contour right after the throat, with little losses as the region feels a strong favorable pressure gradient, permitting a rapid expansion without inducing separation. Bell nozzles areas the most used nozzle design in rocketry due to good performance and reliability results [3].

Minimum length: The most straightforward design for a bell nozzle has a zero-degree inclination at the exit, assuming a philosophy in which 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 7, is defined by a constant 𝜃=0º and 𝜈 =𝜈(𝑀_𝑒). Node 34 represents the intersection of the last right characteristic with 𝐾-=𝜈(𝑀_𝑒), originating from the lip at A, where, owing to the centered expansion fan, 𝜃=𝜈. Utilizing Equation 16, this results in 𝜃_𝑎=𝜃_{𝑀 𝐴𝑋}=0.5𝜈(𝑀_𝑒). 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º<𝜃≤𝜃_{𝑀 𝐴𝑋}, until they intersect the nozzle’s axis and undergo reflection. Subsequently, the points on the last right characteristic line are projected until they intersect 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 [8]. Mishra [19] 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 arises from the MOC’s omission of isentropic losses and the BL.

Corrections applied to the minimum length contour: 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 mandatory, with the calorically perfect assumption inducing significant errors. The MOC can be implemented, but with a local vibrational equilibrium assumption, with a thermically perfect gas model, where the adiabatic index is a function of temperature as per Equation 40, with Θ=3055.556𝐾 [20]. The improved method can provide a more reliable contour for both rocket nozzles and supersonic wind tunnels [21]. However, for hypersonic wind-tunnels, other approaches are needed, as for 𝑀>8 dissociation becomes relevant, with the thermically perfect gas delivering inaccurate results. The in-house algorithm was modified to consider 𝛾 as a function of temperature. An iterative method was implemented, utilizing the Mach number of the respective node.

With the process initiating with 𝛾₀=𝛾_{𝑝𝑒𝑟𝑓}=1.4, it assumes isentropic expansion from stagnation with 𝛾=𝛾_{𝑖𝑡𝑒𝑟𝑎𝑡𝑖𝑜𝑛}, providing a temperature that is used to calculate 𝛾 for the next iteration. Special attention needs to be given to the Prandtl-Meyer angle inversion [22], as it also becomes part of the iterative process. Since it’s assumed a thermally perfect gas, temperature is no longer given as an isentropic expansion. Larger values of stagnation temperature led to longer nozzles, as a higher temperature represents lower 𝛾 values, which, from calorically gas studies, showed that higher area ratios and longer nozzles are needed. There are some exceptions, as there’s a shorter nozzle for 𝑀_exit=4.5 at a stagnation temperature of 1500𝐾. For the same stagnation temperature, a bigger Mach number at the exit has a larger 𝛾 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𝐾 for 𝑀_exit=4.5 with 𝛾_exit=1.4, the dimensions of the contour are slightly different from the perfect gas one, showing some influence of the thermally perfect gas in the modeling of the kernel. Sun [23] 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.

Figure 7:Minimum length bell nozzle design [8].

The traditional MOC also neglects the crucial influence of the BL. This thin region near the nozzle wall, where viscous forces are significant, impacts the flow’s 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 BL. The BL’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 [24]. The presence of a BL introduces a displacement thickness, as shown in Equation 41. 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. Furthermore, the flow field and BL exhibit a two-way interaction, impacting the final exit velocity profile. Higher velocities and stagnation temperatures lead to a thicker BL, making its influence non-negligible in high-enthalpy flows. Consequently, the inviscid contour requires displacement, as detailed in Equation 42 [25].

There are simple algorithms available to predict the BL, and commercial CFD tools can be utilized for iterating the final contour. Two essential parameters characterizing BLs are the shape factor 𝐻, which indicates whether the BL is expected to be laminar or turbulent, and the momentum thickness 𝜑, defined in relation to the momentum flow rate within the BL. These three parameters are interrelated as described by Equation 43. The accuracy of modeling a BL relies on understanding the velocity profile within the BL. Equation 44 presents the von Kármán momentum integral equation for planar compressible flow. Various simplified methods can be employed to solve the momentum differential equation. For instance, applying a Stewartson transformation to 𝐻 and 𝜑 and substituting them into Equation 44. 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 𝐶_𝑓 and incompressible equivalents.

Another approach involves using turbulent BL 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 45, and the BL displacement thickness is determined by Equation 46. The Reynolds number is defined by Equation 47, and viscosity is calculated using Sutherland’s formula [26] for nitrogen gas, as shown in Equation 48.

Ding [27] employed both methods to simulate the BL 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. Still, the simplification of considering turbulent BL over a flat plate under-predicts the displacement by around 35%. Although the BL theory for a flat plane is less accurate, due to its simplicity it was the method chosen, as it provides an understanding of the BL. An in-house algorithm was implemented for BL correction to the calculated contour assuming both thermically and calorically perfect gas. Figure 8 illustrates the inviscid and BL corrected contours, both for thermically perfect gas, with parameters set to 𝑀_exit=2.4, 𝑇_𝑡=1500𝐾, 𝑃_𝑡=1𝑀𝑃𝑎 and 200 characteristic lines. The analysis revealed that, for all cases, the thermically perfect gas assumption resulted in a thicker BL due to the lower adiabatic indexes allowing for lower Reynolds numbers, as per Equation 47. An increase in the stagnation temperature of the main flow was found to increase the BL thickness, while decreasing the stagnation pressure produced a similar effect. The application of the BL displacement thickness method models the development of the BL, and 𝛿∗ at the exit is larger for longer nozzles, even with a higher Reynolds number. A simplification of this method involves neglecting the BL at the throat, which can be substantial and crucial for designing a nozzle with the proper area ratio.

Figure 8:Minimum length bell nozzle with bl correction for thermically perfect gas.

Parabolic contour: Rao [28] 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. A parabolic nozzle consists of a bell-shaped contour where the initial expansion section is defined by a circular arc, followed by a straightening section traced as a parabola. 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 𝜃_𝑁 and the exit divergence angle 𝜃_𝐸. Additionally, the equivalent cone nozzle, and the respective fraction f can be inputs. An in-house code was developed to draw a parabolic bell nozzle as a Bézier quadratic curve [29].

It begins by defining the exit Mach number, determining the area ratio as in Equation 6, and specifying the reference cone nozzle of half-angle 𝛼 and length fraction f. The length and height of the parabolic nozzle are then estimated. Once, 𝜃_𝑁 and 𝜃_𝐸 are stipulated, the parabola can be traced. To predict the performance of the precited contour, another in-house algorithm initiates the MOC at a vertical line, set at the abscissa of the sonic point of the axisymmetric 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. Figure 9 shows the typical net of characteristics and nodes used to evaluate the parabolic nozzle performance at the design point. 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. Since it is a weak shock, theoretically exhibiting almost isentropic behavior, no more sophisticated analyses than the MOC were be employed. 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 cone 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 [29]. 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 particular ly 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 [15].

Figure 9:Parabolic bell nozzle for 𝑀_𝑒=2.4, 𝛾=7/5, 𝜃_𝑁=30°, 𝜃_𝐸=15° and 80% of a 15° cone nozzle.

Truncated ideal compressed contour: The TIC nozzle is derived from an ideal contour, traced in a similar method of the minimum length nozzle, but with a circular arc throat and with the contour traced with mass conservation considerations. It’s important to note the design Mach number of the ideal nozzle must be greater than that of the desired for the TIC nozzle. The ideal contour is then truncated at the abscissa with the ordinate equal to the design exit area ratio for the TIC nozzle. Further, a compression factor is linearly applied to all points of the truncated contour. The previous algorithm is illustrated by Figure 10. The compression introduces a discontinuity in the slope of the curve at point A, as the angle 𝜃 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, with possible changes to the area ratio of the nozzle. Further consequences of the correction are more aggressive initial expansion, 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. Studies reveal 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 contours 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 [30]. An in-house code was developed to trace TIC contours, posing the disadvantage that designs are only possible 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 exi

Bell nozzles for wind tunnel applications

Wind tunnels are crucial elements in the aerospace industry and research, as they offer relatively low-cost means of simulating real events experienced by vehicles. With the advent of the space exploration era and advanced 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 parallel uniform flow. In wind tunnel applications, it is crucial that none of the sections have discontinuities, as this would lead to divergences and velocity non-uniformity at the test section. In the 1930s and 1940s, attempts to trace streamlines through geometry methods were introduced, but they proved to be inefficient and often introduced inaccuracies. A faster and more practical analytical method was later presented by Foelsch [31]. Sivells created a sophisticated method formalized in the algorithm CONTUR, already incorporating a BL correction, also based on the MOC. For hypersonic tunnels, were BL correction and varying 𝛾 are insufficient to model high-enthalpy flows, the divergent is often designed using Sivells’ method with further optimization through commercial CFD tools, preferentially involving chemical dissociation modulation. Potter [32] reported the design of two wind tunnels for testing 𝑀=9 and 𝑀=10, utilizing Sivells’ method with some modifications based on the method proposed by Cresci and employing a BL correction algorithm. However, contrary to expectations, the recorded Mach numbers at the test section were 𝑀=10.15 and 𝑀=9.3, respectively. Benton [33] discovered through experimental testing that for high Mach numbers, BL correction algorithms fail to simulate BL thickness with lengths close to the transversal section diameter. For 𝑀=13.5 and 𝑀=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 a Euler Equations solver, disregarding viscosity, allowing for a considerable test section diameter with uniform and parallel flow [34].

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.

Figure 10:TIC nozzle design steps.

Foelsch method: In a time when only graphical methods were available, Foelsch [31] proposed an 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 requiring some equations derived from it. Figure 11 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 𝜃_𝑃, 𝑟_𝑝, and 𝑀_𝐴≤𝑀_𝑃≤𝑀_𝐸. 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 𝑀_𝑃, the corresponding point K on the contour is given by Equations 49 and 50.

Figure 11:Contour to convert radial flow into parallel flow [31].

The initial expansion zone can be approximated as an arc-circle, as seen in Figure 12, with radius 𝑅, followed by a straight line, with length 𝜆, until meeting the desired area ratio for point A. To ensure a more user-friendly algorithm, 𝑅=𝜆 obtaining Equation 52 and the cone angle 𝜔 will be the slope of the straight line. The circular arc is traced with Equation 53 varying 0º≤ 𝑥 ≤ 𝑅 sin(𝜔). 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 𝑥0, defined in Equation 54.

Figure 12:Contour to convert radial flow into parallel flow [31].

An example of a Foelsch nozzle can be seen in Figure 13. 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 14. A smoother initial curve can be obtained by changing the position of the beginning of the transition curve, as per Equation 55, and from there, tracing a new initial curve, as in Equation 56. The results of this corrected algorithm are seen for the same nozzle in Figure 15.

Figure 13:Foelsch nozzle from in-house code for 𝑀_𝐸=2.4, 𝛾=7/5 and 𝐷=3.15.

Figure 14:Foelsch nozzle from in-house code for 𝑀_𝐸=2.4, 𝛾=7/5 and 𝐷=5.

Figure 15:Foelsch nozzle from in-house code for 𝑀_𝐸=2.4, 𝛾=7/5 and 𝐷=5 with smooth initial curve.

Sivell’s method: 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 [35] 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 16, a transonic solution is established at the circular throat, in this case, the one reported by Hall [5]. 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 𝜔 with the axis. The characteristic lines GE and AB respect the assumption of Equations 57 and 58. 𝜔 values should not exceed 15º and are typically around 12º. Lastly, from B to C, the axial Mach number respects a fifth-degree polynomial, with CD being a straight line with 𝜈=𝜈(𝑀_{𝑒𝑥𝑖𝑡}) and 𝜃=0º. 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.

Figure 16:Sivell’s nozzle scheme [35].

Instead of building an in-house algorithm from scratch, Sivells’ Fortran code, named Contur [36], was utilized. A Python interface, as a library [37], was used to connect the in-house algorithms for further use of this code. There’s an inherent limitation in the inputs as only contours for 𝛾=7/5 can be traced, which 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 17 shows a Sivells contour produced in Contur, with and without BL correction. At first glance, it appears lengthier than the equivalent Foelsch nozzle but with a smoother contour.

Figure 17:Inviscid and BL corrected Sivell’s nozzle contour from CONTUR with 𝑀_{𝑒𝑥𝑖𝑡}=2.4

Dual-bell nozzle

First proposed in 1949 by Cowles and Foster [38], dual-bell nozzles are a very interesting concept that delivers altitude adaptation behavior by varying the effective expansion ratio of the flow, discretely. The 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 [39], comparable to that of an Ariane 5 equipped with a dual bell nozzle, a 72% payload increment, from an original 2000kg, can be obtained. Ferrero [40] finds a 1500kg payload increment, already taking into account modern missions and contemporary understanding of the dual bell nozzle’s flow behavior. Taylor [41] 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 [42]:
A. Mode 1 (low altitude operation): Due to high wall pressure at the extension, aligned with the overextended flow from the inner nozzle, there’s 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.
B. Mode 2 (high altitude operation): Lower wall pressures at the extension and flows under expansion make the flow 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 but some additional losses are found. 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. 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. The effect of this aspiration drags 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 [42]. Recent studies have focused on both delaying and accelerating the transition between modes to minimize the impact of side-loads [39]. 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 [43]:
a) Favorable ΔP<0: Similar to what’s observed in traditional bell nozzles, decreasing pressures will lead to uncontrolled attachment and reattachment.
b) Null ΔP=0: In theory, it promotes an instantaneous transition as soon as the wall pressure, as a whole, drops below the critical value.
c) Adverse Δ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 [41] 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 [43] 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 [40] 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 [44] 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 drags 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 under expansion, near the exit lip, some zones with favorable pressure gradients were observed, which could jeopardize a perfect transition. Kbab [45], [46] and Wu [47] 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 BL correction to the contour provided more accurate results [45].

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 [41] and Kbab [45] 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 [48], but with an arc-circle downstream of the geometrical throat, with 𝜃_𝑁1=30º, as it is a typical value for traditional bell nozzles. Depicted in Figure 18, it consists of two top nozzles. Both nozzles are referenced to a cone nozzle with a half angle of 𝛼=15º, a length corresponding to 𝑓=0.8, and accepting an exit divergence angle of 𝜃_𝐸1=𝜃_𝐸2=15º. 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 59.

Figure 18:Dual Bell nozzle contour from in-house code.

Otsu [48] investigated the acceptable values of 𝛽, 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 𝛽≤1 led to long transitions between modes, starting at higher ambient pressures. For 𝛽=1.2 and 𝛽=1.5, not only was the transition fast, taking only 3ms, but it also occurred at a lower ambient pressure, closer to the ideal value. The depicted nozzle has 𝛽=1.2 and 𝜃_𝑁2=26.05º, representing a good trade-off between a swift transition and manufacturing feasibility.

Aerospike

In the effort to achieve an altitude compensating nozzle, the aerospike has been proposed as an “inside-out” bell nozzle, where the exhaust gases are expanded with a free boundary to the ambient, allowing the adaptation to the changing pressure without constrains nor moving parts [49]. At the cost of some thrust, the spike can be truncated in order to have a more weight efficient nozzle and avoids some extreme heat fluxes [42]. Special attention must be given to the fact that besides the thrust produced by Equation 1, the force that the gases exert on the spike needs to be considered, being represented by Equation 60. For truncated spikes, the base area pressure profile also results in force as peer Equation 61 [50].

H Greer [51] developed a technique for quickly designing an aerospike nozzle contour by utilizing the centered expansion fan and geometric relationships that connect streamlines with solid surfaces. This approach involves plotting a given streamline to serve as the spike’s contour. The initial assumption for this method is that the flow at the throat is sonic and follows a straight line. As illustrated in Figure 19, another assumption is that the flow expansion is solely due to the centered expansion fan at lip A. Mach lines expand and align the flow to achieve the desired exit Mach number with full parallelism to the plug centerline. Lastly, to ensure the flow is straight at the exit, it is assumed that the flow at the throat has an incident angle equal to the Prandtl-Meyer angle of the selected exit Mach number. G Angelino [52] expanded upon Greer’s work by developing a simpler method that calculates contour points discretely, enhancing accuracy. This method involves calculating the length of the Mach line from lip A to the streamline passing through lip B for a range of Mach numbers from 1 to the desired exit Mach number, as well as determining the inclination angle relative to the horizontal. Angelino’s approach involves two steps for each Mach number, to find the position of the corresponding contour point: determining the distance from lip A to the contour streamline and its respective inclination. The length of this line for a given Mach number can be obtained as a dimensionless value based on the throat length, using Equation 62, which relates the throat area to the area crossed by the velocity timeline for that Mach number.

Figure 19:Sketch of the theoretical flow field of an aerospike..

This can be further simplified to Equation 63, with 𝜖(𝑀) defined in Equation 6 and 𝜆 representing the dimensionless length. The inclination is given in Equation 64, derived from geometric relations. Angelino’s simplified method yields a nozzle contour very similar to that obtained from the MOC, with only slight deviations at the nozzle’s end, which can often be truncated. The pressure ratio between the flow and the stagnation pressure is also consistent, with larger deviations occurring at the spike’s end, which can also be truncated in many designs. Although the simplified method involves calculating discrete points to be connected, more complex methods, such as treating points immediately after the throat as part of a circumference, do not produce significantly different contours compared to simple interpolation of all points. Wang [53] demonstrates that Angelino’s simplified method aligns well with results from more complex commercial numerical analysis tools and hot fire test experiments.

E-D nozzle

Rao [54] proposed the concept of the Expansion-Deflection (E-D) nozzle, consisting of an annular asymmetrical throat that is tilted towards a contoured wall. The gases are expanded around a central plug or pintle, and the contour ensures they maintain a mostly axial direction at the exit plane. 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, but with the expansion controlled from inside the nozzle, as the throat delivers flow outwardly, meaning there’s an objective area ratio ceiling. Figure 20 shows the flows field within 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. 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. Losses 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 [42]. 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 [55] 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 [41]. The high viscosity effects felt during open wake are confirmed by Schomberg [56], 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 [57] 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 Figure 20 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 [57].

Figure 20:Flow field of an E-D nozzle at: a) sea-level operation (open wake) and b) vacuum operation (close wake) [42].

MGN

Overall performance can be optimized with the use of a Multiple Grid Nozzle (MNG), when designed with multidisciplinary methods that allow a lighter solution. By activating or shutting down “nozzlettes”, the effective area ratio can be controlled and so the thrust produced. Controlling the area ratio also brings altitude adaptative capabilities. A disadvantage of the MNG is the high ablative erosion, specially to the throat, which will drastically decrease performance during operation [58]. Each “nozzlettes” can be design utilizing some of the previous described methods for bell nozzle contours. 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.

The review performed on the contemporaneity of the MOC in the design methods of supersonic nozzles and performance prediction, was accompanied by the implementation of several algorithms, all described throughout the manuscript. Further exploration of such algorithms corroborated finds such as cone nozzles never operate in their design, BL correction and thermically perfect gas affect the contour, especially for high enthalpy flows and shocks that appear in parabola bell nozzles. Although some more of advanced computational techniques are referred, the MOC offers a unique advantage by transforming complex PDEs into more manageable ODEs along characteristic lines, thus facilitating detailed analysis of supersonic flow fields. The MOC proves a friendly and reliable toll for preliminary assessments where computational resources are limited. Most case studies in literature of nozzle design begin by applying the MOC to draw the first contour. By continuing to refine these methodologies, the aerospace industry can achieve greater performance and reliability in rocket propulsion systems.

This research was funded by the Portuguese Foundation for Science and Technology, I.P. (FCT, I.P.) FCT/MCTES through national funds (PIDDAC), under the R&D Unit C-MAST/Center for Mechanical and Aerospace Science and Technologies, reference: Projects UIDB/00151/2020.

1. El-Sayed AF (2016) Rocket propulsion. Fundamentals of Aircraft and Rocket Propulsion. Springer, London, UK, pp. 907-991.
2. Sutton G, Biblarz O (2001) Rocket propulsion elements. (9^th edn), Wiley, USA.
3. Khare S, Saha UK (2021) Rocket nozzles: 75 years of research and development. Sādhanā, 46(2).
4. Hagemann G, Immich H, Van Nguyen T, Dumnov G (1998) Advanced rocket nozzles. Journal of Propulsion and Power 14(5): 620-633.
5. Hall IM (1962) Transonic flow in two-dimensional and axially-symmetric nozzles. The Quarterly Journal of Mechanics and Applied Mathematics 15(4): 487-508.
6. Kliegel JR, Levine JN (1969) Transonic flow in small throat radius of curvature nozzles. AIAA Journal 7(7): 1375-1378.
7. Dutton JC, Addy AL (1981) Transonic flow in the throat region of axisymmetric nozzles. AIAA Journal 19(6): 801-804.
8. Anderson J (1982) Modern compressible flow, with historical perspective. McGraw-Hill Series in Aeronautical and Aerospace Engineering, McGraw-Hill, New York, USA.
9. (2020) Method of characteristics-internal compressible flows-lesson 6. Ansys, Pennsylvania, USA.
10. Zebbiche T (2019) Stagnation pressure effect on the supersonic minimum length nozzle design. The Aeronautical Journal 123(1265): 1013-1031.
11. Restrepo JC, Bolaños Acosta AF, Simões Moreira JR (2022) Short nozzles design for real gas supersonic flow using the method of characteristics. Applied Thermal Engineering 207: 118063.
12. Kola J, Dvo V (2013) Verification of K𝜔 SST turbulence model for supersonic internal flows.
13. Jia R, Jiang Z, Zhang W (2015) Numerical analysis of flow separation and side loads of a conical nozzle during staging. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering 230(5): 845-855.
14. Zhang Y, Chen H, Zhang M, Zhang M, Li Z, et al. (2015) Performance prediction of conical nozzle using Navier-stokes computation. Journal of Propulsion and Power 31(1): 192-203.
15. Lee C, Choi K, Kim C, Han S (2020) Computational investigation of flow separation in a thrust-optimized parabolic nozzle during high-altitude testing. Computers amp; Fluids 197: 104363.
16. Stlund JO, Muhammad KB (2005) Supersonic flow separation with application to rocket engine nozzles. Applied Mechanics Reviews 58(3): 143-177.
17. Munday D, Gutmark E, Liu J, Kailasanath K (2011) Flow structure and acoustics of supersonic jets from conical convergent- divergent nozzles. Physics of Fluids 23(11).
18. Dakwell HM, Badham H (1963) Shock formation in conical nozzles. AIAA Journal 1(8): 1932- 1934.
19. Mishra NK, Kumar A, Ganesh M, Alam MAR (2019) Design of minimum length supersonic nozzle using the method of characteristics. International Journal of Innovative Technology and Exploring Engineering 9(2): 1370-1374.
20. Benson T (2021) Specific heat capacity-calorically imperfect gas.
21. Matsunaga M, Fujio C, Ogawa H, Higa Y, Handa T (2022) Nozzle design optimization for supersonic wind tunnel by using surrogate-assisted evolutionary algorithms. Aerospace Science and Technology 130: 107879.
22. Hall IM (1975) Inversion of the Prandtl-Meyer relation. The Aeronautical Journal 79(777): 417-418.
23. Sun D, Luo T, Feng Q (2019) New contour design method for rocket nozzle of large area ratio. International Journal of Aerospace Engineering 2019: 1-8.
24. Boynton FP (1968) Exhaust plumes from nozzles with wall boundary layers. Journal of Spacecraft and Rockets 5(10): 1143-1147.
25. Wang Y, Jiang Z (2022) Theories and methods for designing hypersonic high-enthalpy flow nozzles. Chinese Journal of Aeronautics 35(1): 318-339.
26. COMSOL Documentation (1998) Sutherland’s Law.
27. Ding F, Liu J, Huang W, Zhou Y, Guo S (2021) Boundary-layer viscous correction method for hypersonic forebody/inlet integration. Acta Astronautica 189: 638-657.
28. Rao GVR (1958) Exhaust nozzle contour for optimum thrust. Journal of Jet Propulsion 28(6): 377-382.
29. Rick Newlands (2017) The thrust optimized parabolic nozzle.
30. Hoffman JD (1987) Design of compressed truncated perfect nozzles. Journal of Propulsion and Power 3(2): 150-156.
31. Foelsch K (1949) The analytical design of an axially symmetric Laval nozzle for a parallel and uniform jet. Journal of the Aeronautical Sciences 16(3): 161-166.
32. Potter JL, Carden WH (1968) Design of axisymmetric contoured nozzles for laminar hypersonic flow. Journal of Spacecraft and Rockets 5(9): 1095-1100.
33. James R Betnton (1989) Design and Navier-Stokes analysis of hypersonic wind tunnel nozzles.
34. Gaffney RL (2006) Design of a pulse-facility nozzle using the rotational method of characteristics. Journal of Spacecraft and Rockets 43(6): 1216-1224.
35. Sivells JC (1970) Aerodynamic design of axisymmetric hypersonic wind-tunnel nozzles. Journal of Spacecraft and Rockets 7(11): 1292-1299.
36. Sivells JC (1978) A computer program for the aerodynamic design of axisymmetric and planar nozzles for supersonic and hypersonic wind tunnels, DTIC, p. 151.
37. Noah Stockwell (2022) ConturPy: Python hypersonic nozzle design.
38. Cowles FB, Foster CR (1949) Experimental study of gas-flow separation in overexpanded exhaust nozzles for rocket motors.
39. Immich H, Caporicci M, Immich H, Caporicci M (1997) Status of the FESTIP rocket propulsion technology programme. 33^rd Joint Propulsion Conference and Exhibit, American Institute of Aeronautics and Astronautics, USA.
40. Ferrero A, Conte A, Martelli E, Nasuti F, Pastrone D (2022) Dual-bell nozzle with fluidic control of transition for space launchers. Acta Astronautica 193: 130-137.
41. Taylor N, Steelant J, Bond R (2011) Experimental comparison of dual bell and expansion deflection nozzles. 47^th AIAA/ASME/SAE/ASEE Joint Propulsion Conference Exhibit, American Institute of Aeronautics and Astronautics, USA.
42. Hagemann G, Immich H, Van Nguyen T, Dumnov GE (1998) Advanced rocket nozzles. Journal of Propulsion and Power 14(5): 620-634.
43. Frey M, Hagemann G (1999) Critical assessment of dual-bell nozzles. Journal of Propulsion and Power 15(1): 137-143.
44. Wang Y, Lin Y, Eri Q, Kong B (2022) Flow and thrust characteristics of an expansion-deflection dual-bell nozzle. Aerospace Science and Technology 123: 107464.
45. Kbab H, Sellam M, Hamitouche T, Bergheul S, Lagab L (2017) Design and performance evaluation of a dual bell nozzle. Acta Astronautica 130: 52-59.
46. Kbab H, Abada O, Haif S (2023) Numerical investigation of supersonic flows on innovative nozzles (Dual Bell Nozzle). Journal of Applied Fluid Mechanics 16(4): 819-829.
47. Wu K, Sohn GC, Deng R, Jia H, Kim HD, et al. (2023) Study on wall pressure and hysteresis behaviors of a novel dual-bell nozzle. Journal of Mechanical Science and Technology 37(9): 4639-4646.
48. Otsu H, Miyazawa M, Nagata Y (2005) Design criterion of the dual-bell nozzle contour. 56^th International Astronautical Congress of the International Astronautical Federation, the International Academy of Astronautics and the International Institute of Space Law, American Institute of Aeronautics and Astronautics, USA.
49. Nair PP, Suryan A, Kim HD (2017) Computational study of performance characteristics for truncated conical aerospike nozzles. Journal of Thermal Science 26(6): 483-489.
50. Scott J (1999) Aerospike aerodynamics.
51. Greer H (1961) Rapid method for plug nozzle design. ARS Journal 31(4): 560-561.
52. Angelino G (1964) Approximate method for plug nozzle design. AIAA Journal 2(10): 1834-1835.
53. Wang CH, Liu Y, Qin LZ (2009) Aerospike nozzle contour design and its performance validation. Acta Astronautica 64(11-12): 1264-1275.
54. Rao GVR (1960) Liquid rockets and propellants. American Institute of Aeronautics and Astronautics, USA.
55. Taylor NV, Hempsell CM (2004) Optimizing expansion deflection nozzles for vacuum thrust. The Aeronautical Journal 108(1088): 515-522.
56. Schomberg K, Doig G, Olsen J (2014) Computational simulation of an altitude adaptive nozzle concept. Applied Mechanics and Materials 553: 223-228.
57. Paul PJ, Nair PP, Suryan A, Martin MJP, Dong Kim H (2020) Numerical simulation on optimization of pintle base shape in planar expansion-deflection nozzles. Journal of Spacecraft and Rockets 57(3): 539-548.
58. Chasman D, Birch M, Haight S, Graffam R (2005) A multi-disciplinary optimization method for Multi Nozzle Grid (MNG) design-final report. 43^rd AIAA Aerospace Sciences Meeting and Exhibit, American Institute of Aeronautics and Astronautics, USA.