A shape design optimization methodology based on the method of characteristics for rocket nozzles

Even though shape optimization is a powerful tool for designing aerospace vehicles, it can be time-consuming when high-fidelity models are employed. Thus, lower-fidelity simulations covering a wider design space can be a solution for shape optimization in the early design phases. With this in mind, the present work aims to develop a low-fidelity and fast method to conduct nozzle shape optimization. This method consists in using the free-form deformation (FFD) parameterization technique to control the nozzle shape by means of an optimization algorithm to maximize the coefficient of thrust determined by a two-dimensional method of characteristics (MoC). To verify the reliability of the proposed method, a similar optimization process is carried out, recurring to high-fidelity simulations, namely using an Euler solver, in the open-source framework SU2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {SU}^2$$\end{document}. This latter optimization process is established as a surrogate-based optimization (SBO) not only to mitigate the SU2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {SU}^2$$\end{document} framework limitations in performing shape optimization on nozzles, but also as a way to reduce the computational power. A good agreement between the results from both methods is achieved, displaying solely a small offset concerning the optimal contour width and the coefficient of thrust. Hence, this proves the usefulness of the developed shape optimization strategy based on the MoC for the preliminary design of nozzles.


Introduction
The universe and the vastness of space have always been one of the greatest mysteries known to mankind.As time passed and technology evolved, new secrets about the cosmos were revealed and humanity's interest in space grew exponentially.
As of September 2021, more than 4500 satellites orbit the Earth simultaneously [1].This number is only expected to grow higher and higher without any predictions to stop, since most communication technologies of this day and age rely on satellites.Even though the operational cost of rocket launchers has been declining in recent years due to improvements in their re-usability, space flight is still only available to governments, major companies and hyper-wealthy persons.
Having a great impact on the overall performance of a rocket launcher, the proper assembly of a rocket nozzle weighs immensely on the amount of fuel needed for a specific mission [2], and consequently, the amount of payload supported, since the saved propellant weight can be replaced by more payload [3].Therefore, the more optimal a rocket nozzle is, the less money has to be invested into a space launch.
Nozzles were first introduced to build a device capable of changing a flow's characteristics, such as its velocity and pressure distribution.Carl Gustaf Patrik de Laval developed the first convergent-divergent nozzle capable of increasing a jet flow into a supersonic state in 1890, which later became known as a "de Laval" nozzle [4].Robert Goddard is credited with the first flight using a "de Laval" nozzle with a combustion chamber [5].
Rocket nozzles come in a wide variety of configurations, like an ideal, conical, bell, plug, expansion-deflection (E-D) and dual-bell nozzles.Recently, a multi-nozzle grid was also developed.Figure 1 provides illustrative examples of different nozzle designs.
An ideal nozzle is characterized by producing an isentropic flow (i.e., without internal shocks) and having uniform properties at the exit, meaning constant pressure, temperature and velocity over the whole exit plane.However, the ideal nozzle that gives maximum thrust performance is heavy and lengthy.To expand a flow to match near vacuum conditions, one would need a near infinite nozzle [6].The length of an ideal nozzle can be decreased by permitting all the expansion to occur just at the throat (amid a sharp corner) and then constructing the nozzle contour to turn the flow such that it can attain an axial uniform flow at the exit [7].This nozzle is referred to as a minimum-length nozzle.
The application of conical nozzles was very popular during the early stages of rocket propulsion due to its simplistic design and ease of manufacture.However, its main downfall compared to an ideal nozzle comes from the fact that the exit flow is divergent.
Bell nozzles are the most commonly used shapes in rocket engines.They comprise a high-angle expansion section ( 20 • to 50 • ) downstream of the nozzle throat, followed by a gradual reversal of the nozzle contour slope so that at the nozzle exit the divergence angle is small, avoiding major divergence losses [3].Rao [8] developed in 1958 a method by using the calculus of variations to design the wall contour of the optimum thrust nozzle by using a simple parabolic approximation.
The main characteristic of a plug nozzle is its capacity to interact with the external ambient, resulting in having a free jet boundary that acts as a virtual outer wall and expands and compresses to match the free-stream ambient pressure.The separation of flow can be avoided, implying that plug nozzles are altitude compensating [9].
Like the previous nozzle configuration, the E-D nozzle is altitude compensating.Although looking similar to the bell nozzle, the flow is turned by a "center body" onto the outer diverging nozzle wall [3].This results in the creation of a viscous wake region within the nozzle.
The dual-bell nozzle, as mentioned in the name, consists of two distinct contours between the throat and the exit [10].Considered an altitude-adaptive nozzle concept, it combines a small nozzle area ratio at low altitude with a large one at high altitude.Its design is typically based on an inner base nozzle, followed by a wall inflexion region and an outer nozzle extension.
The multi-nozzle grid configuration [3] is the most recent and the only nozzle configuration where the nozzle length (i.e., plate thickness) to throat diameter can be less than one, yet is capable of providing an extremely high area ratio.
Traditionally, rocket nozzles have been designed exclusively for propulsive performance.The optimization of nozzle profiles for maximum thrust was pursued long before computational fluid mechanics (CFD) or optimization became widely available tools.Some examples combining these two numeric tools can be found in [11,12] and [13].
Classical optimization procedures usually began with an inviscid design, such as Rao's method [14].Then a boundary-layer correction would be added to compensate for the viscous effects.Recent advances in computational technology have allowed the integration of the full Navier-Stokes equations, as well as enabling automatic design methods developed by combining the CFD and optimization codes [15].
Numerical methods can for the most part be divided into low-and high-fidelity methods.Both are extensively used because of the computational cost of higher-fidelity models.When computationally expensive analyses are used, it is often of benefit to utilize multiple levels of fidelity to reduce the size of the design space and initialize subsequent analyses from a lower-fidelity solution [16].This approach is referred to as the multi-fidelity approach.
In [17] and [18], Asha et al. presented a study about the design of a nozzle with a minimum length using the method of characteristics (MoC).As a conclusion, both articles recommend the MoC as the most appropriate to be used for this supersonic nozzle design.
Matsunaga et al. [19] published a paper presenting a design methodology for a supersonic wind tunnel nozzle with the objective to obtain a uniform flow at the nozzle exit.To archive this objective, they conducted a shape optimization of the nozzle using a surrogate model based on CFD simulations and use an evolutionary algorithm to identify the global optimum solution.Even though this application was performed for a wind tunnel nozzle, it provides interesting aspects for a throttle nozzle design.With this motivation in mind, the main objective of the current paper is to develop a low-fidelity model using a MoC capable of being integrated in an optimization framework to provide accurate and fast results during the preliminary design phase of rocket nozzles.
To accomplish this goal, different features from the aforementioned methods are combined in an optimization framework that constitutes the novelty of the paper.This framework consists in optimizing the nozzle external shape parameterized by means of free-form deformation using the MoC.Its performance is then compared to a similar optimization strategy based on a surrogate model built from CFD simulations.
This article is ordered into five sections.Each one is subsequently divided into subsections aiming for a clear organization and smooth reading.Firstly, Sect. 1 describes the objectives and motivation behind the completion of this research and presents a review on the topics of rocket nozzles and their optimization.The most relevant mathematical formulations related to the field of fluid mechanics, optimization and parameterization are presented in Sect. 2.Then, Sect. 3 employs the mathematical foundations previously reviewed and describes the steps and strategies toward achieving the results.This is followed by Sect.4, where the results regarding both low-and high-fidelity-based optimizations are displayed.Finally, the main conclusions regarding the obtained results are drawn in Sect. 5.

Nozzle flows
Nozzle flows describe all types of flows going through nozzles, which are devices that increase the velocity of a fluid at the expense of pressure.Their interior flows experience fascinating phenomena ranging from supersonic speeds to shock waves.
Assuming quasi-one-dimensional flow properties, it is possible to obtain an equation that relates the change in velocity du to the change in area dA assuming isentropic flow [20], where M is the Mach number.It is important to notice that to reach sonic flow, one must first accelerate the flow by decreasing the nozzle area.After achieving sonic conditions, the area A can be increased to further increase the flow's velocity u.Hence, a nozzle designed to achieve supersonic flow at its exit is a convergent-divergent duct.The minimum area, where sonic flow is achieved, is called the throat.
Another very important equation is the area-Mach number relation, which is deduced in [20]: where is the ratio of specific heats.This equation shows that the Mach at any point of the nozzle is a function of the ratio between the local and the throat area.

Prandtl-Meyer expansion
Expansion fans occur when a supersonic flow is turned away from itself [20].An expansion fan is a continuous expansion region that can be visualized as an infinite number of Mach waves, contrary to an oblique shock that is made of a single shock wave, a consequence of the flow turning into itself.An expansion wave appearing from a sharp convex corner is called a centered expansion wave.This kind of fan will be visualized during the implementation of the Method of Characteristics and is commonly addressed to as Prandtl-Meyer expansion waves in honor of Ludwig Prandtl and Theodor Meyer, who first worked on a theory for centered expansion waves.From Prandtl-Meyer's theory comes an equation that relates the infinitesimal change in velocity dV to the infinitesimal deflection d : For a convex corner of angle d Eq. 3 is integrated resulting in: where is called the Prandtl-Meyer function:

Method of characteristics
When dealing with supersonic flows, the partial differential equations (PDEs), i.e., the Euler equations, that govern the flow-field properties (e.g., temperature, pressure, speed and density) are hyperbolic.The main goal of this method is to simplify these PDEs into an ordinary differential equation (ODE).Characteristic lines divide two adjacent regions described by expressions which are analytically different.In general, characteristic lines exist when the transition from one region to another involves discontinuities of some derivatives.This means that there are certain lines (2) A A t in the xy space along which the derivatives of the flow variables are indeterminate.These characteristic lines divide the region of the flow where the action is produced from the region of the flow that ignores the presence of the disturbances [21].
When the flow is governed by a supersonic steady motion having constant entropy and constant total enthalpy, the characteristic surfaces are coincident with the envelope of Mach as shown by the following equation: where and u = Vcos( ) and v = V sin( ).
Along the characteristic lines displayed in Fig. 2, the PDE reduce to an ODE, as presented next: Equation 8 is identical to the expression used to describe Prandtl-Meyer expansion waves resulting in the following relations: where C + and C − denote the left-and right-running charac- teristics, respectively.

Free-form deformation
The idea behind free-form deformation (FFD) is to enclose a geometry into a parallelepiped of control points that define the parametric domain.A physical analogy for FFD would be to imagine a parallelepiped of clear, flexible plastic undergoing deformation [22].FFD can be based on different principles like Sederber's scheme based on Bernstein's polynomials or NURBS FFD which is more complex.
To apply Sederber's technique, every point of the base geometry has to be converted into the coordinate system (x, y, z) established by the parallelepiped region Let P ijk , i = 0, ..., l, j = 0, ..., m, k = 0, ..., n , represent the movement of the control points from their latticed position, then one computes the new points X new by where

Design optimization
Design optimization deals with the improvement of a structure's properties by improving its shape.In the specific case of a rocket nozzle, properties such as the coefficient of thrust and the specific impulse can be analyzed to improve its performance.Optimization is often confused with the term improvement.Mathematically speaking, it means finding the best possible solution by changing control variables, often subjected to constraints [23].
To find the optimal design, it is necessary to define an objective function whose goal is to differentiate better from worse designs.Properties like the coefficient of thrust and the specific impulse can be quantified using an objective function, where different design variables result in different outcomes.Objective functions can be implemented in various ways.In this work, both the method of characteristics and CFD are used to compute the objective functions of interest.
Most engineering applications are constrained problems.Constraints are used to restrict possible solutions to a feasible region.The constraints can limit the design values directly, denominated bounds, or indirectly limit the results (10) Fig. 2 Left-and right-running characteristic lines through point A, adapted from [20] through equality and/or inequality constraints [23].A typical constrained optimization problem can be formulated as: where x is the vector of design variables, f is the objective function, and lb, ub are the lower and upper boundaries of the design variables, respectively.A, b and A eq , b eq are linear inequality and equality constraints, respectively.

Optimization algorithms
The process of optimization consists in iteratively guessing the optimal value of the design variables until a solution is obtained.No single optimization algorithm is effective for all optimization problems, since each different algorithm has a different strategy to move between consecutive iterations [23].A common classification to characterize optimization algorithms is by dividing them into gradient-free and gradient-based algorithms.
The gradient-free optimization approach is solely based on moving from one point to another if the value of the objective function decreases.A gradient-free algorithm is easier to set up because no additional coding is needed other than the objective function and its constraints.On the other hand, gradient-free methods require a lot more objective function evaluations than gradient-based method when the number of design variables increases.
Gradient-based algorithms take advantage of both the objective function and its derivative with respect to the design variables to converge to the optimum more efficiently.Gradients are used to ensure that the optimizer converges to the mathematical optimality condition.Though they usually need fewer iterations, gradient-based algorithms require the objective function to be sufficiently smooth.In practice, however, they can tolerate discontinuities as long as they are not near the optimum [23].

Surrogate-based optimization
A surrogate model consists in an approximation of a functional output that represents a curve fit to some data.The main purpose of a surrogate-based optimization is to perform optimization using a model that is much faster to compute than the original function, without losing substantial accuracy [23].Surrogate models are a computationally cheap and easy alternative for models which are computationally expensive, such as high-fidelity CFD simulations.
In addition, they are also helpful to visualize how the objective function varies with respect to the design variables.(13) On the contrary, surrogate models display poor scalability; in other words, the larger the number of inputs, the more are the model evaluations needed to construct a surrogate model that is accurate enough.

Implementation
The development of the implementation process presented in the next section can be visualized through a flowchart presented in Fig. 3.

Minimum-length nozzle: method of characteristic implementation in MATLAB®
Supersonic nozzles can be divided into two different types: gradual-expansion nozzles, and minimum-length nozzles.
Gradual-expansion nozzles are typically used when maintaining a high-quality flow at the exit is desired, like for supersonic wind tunnels.However, due to being lengthy and heavy, a minimum-length nozzle, which utilizes a sharp corner at the throat, is a better option for a rocket nozzle [24].
If the nozzle contour is not properly shaped, shock waves can occur inside the duct.The method of characteristics provides a technique for properly designing the contour of a supersonic nozzle for shock-free, isentropic flow, taking into account the multidimensional flow inside the duct.
As mentioned previously, for a minimum-length nozzle, the expansion section is shrunk to a point and the expansion takes place through a centered Prandtl-Meyer wave emanating from a sharp-corner throat with an angle w max as demonstrated in Fig. 4. The criterion for a minimum-length nozzle is determined as half of the Prandtl-Meyer function for the exit Mach number: Another important aspect regarding the minimum-length nozzle is that its contour is made to absorb expansion waves instead of reflecting them.This means that the flow near the wall will neither expand nor compress, meaning its direction angle stays equal according to Eq. 8.The condition for no compressibility is

Minimum-length nozzle: algorithm implementation
The method of characteristics needs some initial data to calculate the flow properties.From the desired exit Mach number M e , using Eq.14 the angle at the starting point is ( 14) obtained.Since all characteristics come up from the same point and have to be initiated with a different initiation angle ini i to display different paths, the initiation angle vector ini is computed as a distribution with the intervals [0, w max ].

Arbitrary nozzle-algorithm implementation
Previously, the MoC was used to define a contour.However, the main goal is to apply the MoC on preexisting contours and calculate its flow-field properties to calculate its coefficient of thrust.
The major differentiation is that Eq. 15 no longer is obeyed and the direction of the flow near the wall wall is equal to the contour's slope (Euler boundary condition): Implementing the method of characteristics on a contour is not as straightforward as creating a contour, because the grid generated by the MoC is not identical to the domain of the nozzle.This happens because the contour is continuous in space, while the grid is discrete.In other words, the algorithm in charge of computing the characteristics generates each left-running characteristic line, one by one, and will only stop when the last characteristic line goes beyond the nozzle contour.This leads to the last left-running characteristic being partially located outside the nozzle.To calculate (16) wall = contour .

Fig. 3 Implementation process fowchart describing all major procedures
Fig. 4 Minimum-length nozzle design, adapted from [18] the flow-field variables at the exit, a "virtual" characteristic is imagined between the last two left-running characteristics, where the flow-field properties can be calculated through linear interpolation.

Coefficient of thrust: method of characteristics
The thrust F of a quasi-one-dimensional nozzle flow is known and calculated as where ṁ stands for mass flow, V e is the exit speed, p e and p amb are the exit and ambient pressures, respectively, and A e denotes nozzle exit area.As mentioned before, this equation is used for quasi-one-dimensional flow, meaning the flowfield variables only vary in x.
To obtain the thrust of the nozzle using the method of characteristic, it is necessary to apply a Riemann sum to the normal component of the exit velocity and pressure p e along the last two characteristic lines and interpolate the values obtaining the average exit pressure and velocity along the "virtual" characteristic.For this problem, a "centered" Riemann sum is applied as: and the average value Proceeding to the interpolation, one obtains: where x wall 1 and x wall 2 are the x-location of the intersection of the characteristics with the contour before and after the exit point, respectively.After the determination of both exit velocity and pressure, the thrust can be obtained using Eq. 17 as well as the thrust coefficient The subscripts 0 and t denote total conditions and throat, respectively.

Design optimization in MATLAB®
Before implementing any kind of optimization method, it is necessary to define an objective function.As mentioned (17) F = ṁV e + p e − p amb A e , (18) before, the main goal of the optimizer is to increase the coefficient of thrust of a rocket nozzle.When using an optimizer for the MoC, one has to take into consideration that the restriction to the optimizer is intrinsic to the method itself instead of being expressed separately.Equation 13 reduces to: Regarding the optimizer, two different optimization methods were applied to the method of characteristics (fmincon and NSGA-II) which are presented next.
fmincon is a gradient-based method from MATLAB ® used to find the minimum of a constrained nonlinear multivariable function identical to Eq. 13.Like most gradientbased methods it is designed to operate on problems where the objective function and its derivative are continuous.This method performs a local search since its characterized by a single starting point and presents a deterministic behavior because it always evaluates the same points given the same initial condition.
The genetic algorithm NSGA-II stands for non-dominated sorting genetic algorithm II [25].Contrary to the previous algorithm, NSGA-II is a gradient-free algorithm that is population based like most genetic algorithms [26].The optimization starts with a set of design points rather than a single starting point, and each optimization iteration updates this set in some way.This algorithm performs a global search, since it does not begin at a specific point and evolves throughout the problem domain.Furthermore, it involves randomness in the population generation.Another important characteristic of this algorithm is associated with the fact that it is a multi-objective optimizer.
As the second objective to be analyzed, the width, or rather, the height of the nozzle is chosen since it shows a direct proportionality to the drag forces subjected to the rocket nozzle.During take-off, minimizing forces opposed to the thrust is critical, with the drag force being one of the most dominant ones, due to the high air density near sea level.
The NSGA-II algorithm is inspired by biological reproduction and evolution using three main steps: selection, crossover, and mutation [23].

Computational fluid dynamics: SU 2 Solver
Computational fluid dynamics or CFD is a very powerful tool for the analysis of systems involving fluid flow, heat transfer and other associated phenomena.From the 1960s onward, CFD techniques have been used in the aerospace industry.CFD leads to substantial time and cost reduction of new designs as well as enables the study of systems where (23) controlled experiments are difficult, such as simulating rocket nozzle supersonic flows [27].
On the other hand, one has to take into consideration its computational cost.When working with a computational budget, it is necessary to be able to find a trade-off between the accuracy and the computational cost of the solver.
SU 2 is an open-source collection of software tools writ- ten in C++ and Python for the analysis of PDEs and PDEconstrained optimization problems [28].
To properly run the simulation applied to the nozzle mesh, a configuration cfg file had to be developed.Some of the decisions regarding the choices used in the configuration are presented next.
For the simulation of a supersonic nozzle flow, a Euler solver was applied, to compare this high-fidelity method with the Method of Characteristics, since both are based on the Euler equation.As a result, no turbulence or viscosity models needed to be implemented.The fluid was assumed as an ideal gas with the ratio of specific heats = 1.4 and a specific gas constant of R = 287 J/(kg ⋅ K), the default value for standard air.Regarding the boundary conditions, the Euler condition was applied to the nozzle wall, while the symmetry condition was enforced on the centerline.As for the inlet (throat) of the nozzle, a supersonic inflow Dirichlet condition was applied and the temperature, static pressure and velocity direction were assumed.Concerning the output of the nozzle one simply enforced the static outlet pressure.To compute the gradients of the flow variables the weighted least-squares numerical method for spatial gradients was implemented.Regarding the Courant-Friedrichs-Lewy (CFL) condition, an adaptive CFL number was implemented which could fluctuate between 0.5 and 100.A multigrid with two levels is used to carry out the simulations faster and help convergence.The convergence criteria used to observe if the solver converges is the root-mean-square (RMS).Finally and most importantly, an ROE scheme was used as the convective numerical method.ROE is a modern high-resolution, shock-capturing scheme.The main idea behind this scheme is to determine the approximate solution by solving a constant coefficient linear system instead of the original nonlinear system [29].This scheme has shown to be particularly effective when dealing with the Euler equation and was implemented for this reason.

Grid convergence study
Before commencing with the optimization of the nozzle design, a grid convergence study was conducted to define the optimal mesh discretization.The assessment of the spatial convergence of a simulation involves performing the simulation on two or more successively finer grids.As the grid is refined and the time step is reduced, the spatial and temporal discretization errors, respectively, should asymptotically approach zero, excluding computer roundoff error.A usual method to examine spatial and temporal convergence is known as Richardson's extrapolation (RE) [30].The Richardson's extrapolation value, which estimates the continuum value (value at zero grid spacing), is expressed as: where and f 1 , f 2 and f 3 represent the finest to the coarsest grid, respectively.The grid refinement ratio r for a structured mesh is defined as: where N s is the number of cells on the mesh.
A structured mesh can easily be implemented for a rocket nozzle flow without compromising the solver's accuracy or demanding too much computational power.

Surrogate-based optimization: implementation
After completing the grid convergence study, the surrogate model may be initiated.Having obtained a database of meshes, running the simulation is the next step to gather the training data to build the surrogate model.Implementing polynomial fitting (in MATLAB ® and EXCEL ® ), the surrogate model is computed and later cross-validated with the test data, which was generated simultaneously with the training data.If the surrogate model presents negligible errors regarding the test data, one can assume that the model simulates the SU 2 framework to a certain extent.Finally, an opti- mizer (fmincon and/or NSGA-II) can be applied to obtain the optimal design and compare it to the design obtained using the method of characteristics.
Although SU 2 provides a feature to perform shape design optimization, the computational budget available for this work would not sustain such complex simulations.Furthermore, to use this design capability in SU 2 , one has to adapt the source code such that the nozzle design with the aforementioned features is enabled (e.g., nozzle outputs and flow post-processing), which is outside the scope of this document. (24)

√
N s fine ∕N s coarse

Results
Before starting with the discussion of the results, it is important to underline that a two-dimensional MoC produces a wedge-shaped nozzle with a rectangular outlet [31].Additionally, the nozzles produced with the twodimensional algorithm are calculated for a unit breadth and unit throat height.

Minimum-length nozzle: results
Being the easiest nozzle design to implement, the minimum-length nozzle is a good candidate to test the implementation of the method of characteristics.Since the algorithm applied to arbitrary nozzles is a variation of the one developed for minimum-length nozzles, it is of utmost importance to assure the accurate functioning of this algorithm.
The implementation of the MoC for a minimum-length nozzle with an exit Mach M e = 3 , with a varying num- ber of characteristics, n char shows how quickly the relative error between the numerical simulation and the analytical equation 2 descends by increasing n char .For ten character- istics ( n char = 10 ), this relative error is already under 1%.
Another important aspect of the minimum-length nozzle is noticing that the characteristics maintain an identical distance between each other outside of the kernel zone.This reaches back to Eq. 15, where it was assumed that the flow experiences neither expansion nor compression, as proven by Fig. 5.The kernel zone is defined as the zone of flow expansion and it is outlined by the characteristics that are initialized at the throat.

Bell nozzle nozzle
To design a bell nozzle contour, the parameterization method FFD is applied to a base geometry, which in turn is generated from a minimum-length nozzle.For all bell nozzle simulations, the contour of a minimum-length nozzle with M e = 3 and n char = 50 was chosen as the base geometry.
Before deforming any kind of geometry, a free-form deformation box has to be built around the base geometry.The starting point as well as the end points of the nozzle's geometry delimit the size of the box.A FFD box can have as many control points as wanted; however, the influence of each control point reduces when increasing the refinement of the grid, as shown by Eq. 11 and portrayed in Fig. 6.
The method of characteristics can now be used for different bell nozzles by changing the deformation parameters x FFD , which are linked to P ijk .The node most upright of the FFD box is the one that will be dislocated in the vertical direction to obtain new geometries.This decision is taken, because that control node is the one that has the most influence on the coefficient of thrust, since the exit velocity and pressure are mainly dependent on the exit area, which is directly controlled by that specific node as illustrated in Fig. 6.

MoC-based optimization in MATLAB®: results
Before carrying out any optimization, some properties of the rocket nozzle have to be clarified.The base nozzle is engineered in a way that its exit pressure equals 1 atmosphere, where 1 atm = 101325 Pa.To achieve this feat, the chamber pressure p 0 has to be exactly 3723300 Pa.The noz- zle is structured this way for it to be optimal for an ambient pressure near the sea level, since the condition of optimality states that p e = p amb .Additionally, a chamber temperature T 0 of 3000 K was chosen.
Various simulations were run using the gradient-based method.As discussed before, the most upright node of the FFD box has the biggest influence on the coefficient of thrust of a rocket nozzle.For that reason, this will be the only node to be deformed, with a lower and upper boundary of − 0.5 and 2, respectively.To have some variety throughout the results, four different FFD grids will be used.
The results, displayed in Table 1, show that for pressures lower than sea level, the nozzle widens, as expected since the nozzle tries to minimize the difference between the exterior and the exit pressures (optimality condition).This outcome can be visualized for a [2× 2] FFD grid in Fig. 7.However, a total optimal condition ( p e = p amb ) cannot be achieved because the flow may not be expanded too extensively.
To achieve the best possible results, one can allow other control points to be deformed as well as shown in Fig. 8.By applying this change the C T can be increased even further.
The optimization achieves an overall increase of the coefficient of thrust in vacuum C T vac from1.56784 to 1.61612.Regarding the NSGA-II method, an identical result was attained, as can be seen in sect.4.5.

Grid convergence study SU 2 : results
A grid convergence study is executed for five different structured meshes.The study showed that the C T converges when the number of grid elements is increased as well as the RE estimated, calculated from the three finest grids.For generating the surrogate model, the data set was computed using a [ 200 × 30 ] mesh to keep the computational cost low.Since the computational time increases exponentially with the increase of grid elements and a sufficiently large data set has to be simulated to build the surrogate model, it is wiser to keep the computational cost on the lower end.Nevertheless, it is worth mentioning that a CFD simulation with the chosen mesh still took 3.4s to run, which is substantially higher than 0.12s of an MoC simulation, using the same computer.

Surrogate-based optimization: results
A surrogate-based optimization approach is followed to allow for a fair comparison between MoC and CFD-based optimizations without coding in the SU 2 environment some required features.Table 2 provides a comparison between the two methods for two ambient pressures.The other ambient pressures are omitted for the sake of brevity and due to the fact that they do not change the conclusions drawn for this comparison.
Comparing the results of the SBO, the MoC seems to overestimate the coefficient of thrust.Regarding the optimal deformation, in vacuum, the MoC significantly overshoots the value of x FFD for coarse FFD grids ( Δx FFD ≈ 0.16), while for an ambient pressure of 0.5 atm finer FFD grids are the ones dealing with overshooting issues ( Δx FFD ≈ 60%).For sea-level conditions, the surrogate model outcomes dictate a small deformation, which can be neglected and assumed zero, meaning the nozzle geometry is developed for an optimal nozzle at take-off.The optimal nozzle evolution for various ambient pressure is shown in Fig. 9.
The discrepancies can be caused by various factors.These can be divided between the limitations of the MoC and the limitations of the surrogate model.

Limitation of the surrogate model versus the method of characteristics
The surrogate model was implemented using a wide array of data points, meaning that the polynomial fitting can lead to errors due to "overfitting".Another limitation is due to the value of the data points themselves being poorly precise because the SU 2 framework only outputs the C T with four decimal places.In some cases, the value of C T is equal for adjacent design variables, resulting in difficulties in accurately fitting the data distribution to a polynomial without having to increase its order, inducing the occurrence of "overfitting".As mentioned before, during the grid convergence study, the sizing of the CFD mesh has an impact on the underestimation of the C T .Finally, the main limitation of the sur- rogate model approach is its inability to work with more than one design variable at a time without the curve fitting process becoming too complex.
Being a low-fidelity method, the MoC has considerably more limitations that any high-fidelity model.The main downside of the MoC is its exit conditions not being coincident with the real nozzle exit.This leads to error when the flowfield is not uniform leaving the last characteristic line, since part of the flow is not simulated.This is particularly damaging when the deformation applied to the base geometry only occurs toward the exit (finer FFD grids), since the expansion of the flow near the center of the nozzle is not computed.Another problem leading to a loss in accuracy is the fact that the characteristic grid becomes coarser along the nozzle, reinforcing the statement that variations in the flow due to geometry changes near the exit are less precisely computed.Additionally, the "overshooting" and intersection of the last characteristic line with a made-up contour, as well as the interpolation and Riemann sum, induce errors, since the distance between the "virtual" exit characteristic and its adjacent ones is not constant.
One might ask if all these inaccuracies could be avoided by increasing the number of characteristic lines used during the implementation of the MoC.The problem with that statement is that for too fine a grid the characteristic line might collapse when experiencing slight compression.This outcome is strictly undesired, since the whole characteristic grid might break down introducing far worse errors than the ones trying to be avoided.

Conclusions
The present work consisted in developing a fast and reliable tool for preliminary nozzle design optimization, as well as testing out the reliability of the proposed low-fidelity optimization method.
A simulation tool based on the method of characteristics in 2D was developed.It demonstrates the capability of contouring ideal nozzles with high accuracy.For a total number of 100 characteristic lines, the nozzle exit height h nozzle dif- fers from the analytical value by a margin of 0.18%.
A FFD-based parameterization technique was implemented to automatize the generation of nozzle geometries to be evaluated and optimized.The first major achievement was attaining an automated procedure to deform a base contour and simulate its flowfield.An ideal nozzle with an exit Mach number of 3 was chosen as the base contour, calibrated for optimum thrust at sea-level conditions.
Having established an optimization process, based on the MoC, capable of resembling an objective function, two contrasting optimization algorithms were employed to maximize the thrust (and minimise drag) produced by the nozzle for various ambient pressures.The resulting thrust-optimized contour (TOC) demonstrated being able to be contoured accurately by a parabola, meaning the optimization process deformed an ideal nozzle into a TOP.An improvement of 3.08% is achieved for the coefficient of thrust in vacuum C T vac by the developed optimization process using multiple design variables and a [ 9 × 2 ] FFD grid.
To validate the results obtained by the implementation of the Method of Characteristics a CFD simulation using the Euler solver of the SU 2 framework was configured.Applying a Richardson extrapolation to the various results obtained for C T SSL the study showed a relative error of 0.297% concerning the MoC for an ideal nozzle with M e = 3 .For the following simulations, a mesh of 200 × 30 elements was appointed due to its low computational time of 3.4 s and its still negligible relative difference of 1.157%.However, this run time corresponds to a 96% increase when compared to a single MoC simulation.
An accurate surrogate model based on the CFD training data was built as a way to simplify the optimization process, since implementing a surrogate model-based optimization is considerably computationally cheaper, than running a CFDbased optimization.
Having implemented an optimization based on a highfidelity method, even though in some sort simplified by the implementation of a surrogate model, the fidelity of the MoC-based optimization was verified.A good agreement regarding the results of both routes was achieved, however with a small offset concerning the design variable x FFD and coefficient of thrust C T for the optimal design.The maximal deviation of the design variable Δx FFD between both meth- ods showed not to exceed a 20% absolute variation.
In conclusion, the present work proved that the MoC is a strong and reliable tool for preliminary design optimization since it reduces the computational cost regarding CFDbased optimization for a small loss in accuracy.The method developed throughout this work is ideal for a multi-fidelity approach.

Table 1
Gradient-based optimization results for various FFD grids and ambient pressuresAmbient Pressure p amb = 0 atm p amb = 0.25 atm p amb = 0.5 atm p amb = 0.75 atm Optimal nozzle geometry for [2 × 2] FFD grid for various ambient pressure using a gradient-based optimization approach Fig.8Optimal nozzle geometry for [9 × 2] FFD grid using multiple optimization variables in vacuum

Table 2
Surrogate-based optimization results for various FFD grids and ambient pressuresAmbient pressure Grid size MoC− x FFD SU 2 − x FFD Δ (%) MoC− C T SU 2 − C T Δ (%) Fig. 9 Optimal nozzle geometry for [2 × 2] FFD grid for various ambient pressure using a surrogate-based optimization approach