Multi-Objective Bayesian Optimization Design of Elliptical Double Serpentine Nozzle

: The use of traditional optimization methods in engineering design problems, specifically in aerodynamic and infrared stealth optimization for engine nozzles, requires a large number of objective function evaluations, therefore introducing a considerable challenge in terms of time constraints. In this paper, this limitation is addressed by using a sample-efficient multi-objective Bayesian optimization that takes Kriging as a surrogate model and Expected Hypervolume Improvement as the infill criterion. Using this approach, the probabilistic model is continuously established and updated, and the approximate Pareto front is obtained at a relatively small computational budget. The objective of this work is to evaluate the applicability of employing a multi-objective Bayesian optimization framework for the aerodynamic-infrared shape optimization of an elliptical double serpentine nozzle at 6 km flight condition, where the objective functions are evaluated by means of high-fidelity computational fluid dynamics and reversed Monte Carlo ray tracing simulations. We achieve good results in both infrared radiation signature reduction and aerodynamic performance improvement with a reasonable number of evaluations, indicating that the proposed method is effective and efficient for tackling the computationally intensive optimization challenges in the aircraft design.


Introduction
With the development of infrared (IR) tracking and detection technology, the infrared stealth characteristics of the new-generation fighter will meet new challenges.As an important component of the engine, the nozzle provides the majority of the thrust, while the high-temperature, high-velocity plume and hot parts inside the nozzle make it the primary infrared radiation within the aircraft [1].In order to enhance the survivability and overall combat proficiency, extensive research [2][3][4] has been conducted on their infrared radiation signature and aerodynamic characteristics, which have naturally become the focus of research.
Traditionally, aircraft nozzle typically feature a circular inlet and circular exit with a straight transition of flow from inlet to exit.However, the growing demands for stealth capabilities in fighter aircraft and the drive to reduce noise in civil aircraft have shifted researchers' attention from the conventional axisymmetric nozzle to the asymmetric nozzle [5].The potential of serpentine nozzles for IR suppression has led to their implementation in a series of stealth bombers and unmanned aerial vehicles (UAVs) [6][7][8].The asymmetric shape of these nozzles allows for better mixing between the hot exhaust gas and the surrounding atmosphere, which can result in a shorter high-temperature potential core length of the exhaust plume [9,10].Additionally, the curved design of the nozzle helps hide the major hot sections of the engine, such as the turbine section and central cone, thereby further contributing to the reduction of the IR signature [11].On average, the IR radiation intensity can be reduced by at least 70% compared to a circular nozzle [12].Typically, for the serpentine nozzle, there is often a significant reduction in the IR signature in the axial direction, while relatively higher signature levels are observed on the sides and below the nozzle when the aspect ratio of the nozzle increases [13].
According to the literature mentioned above, it is known that the serpentine nozzle with a non-circular exit can effectively decrease the infrared radiation emitted by the hightemperature walls and exhaust gases.However, the transition from a circular inlet to a non-circular exit, as well as the axial centerline curvature, can induce secondary flows inside the nozzle [5].These secondary flows can significantly complicate the flow field, potentially leading to detrimental effects on the aerodynamic performance of the nozzle.Rajkumar et al. [14] conducted an experiment to analyze the flow characteristics of a shallow single serpentine curved convergent nozzle under different pressure ratios.The results indicated an initial deceleration of the flow due to the turning of streamlines caused by the nozzle curvature, as well as an acceleration of the core flow due to the geometric transition.Furthermore, the vortices induced by separation from the tailcone and the sidewall vortices have a significant impact on the temperature distribution of the nozzle [15,16].This effect becomes particularly prominent in the instance of double serpentine nozzles, where regions with high temperature, often termed as hot streaks, tend to develop along the nozzle surface.These high-temperature zones will cause material degradation of the nozzle and adversely affect the IR stealth capabilities of the nozzle [17,18].
Generally, simulations of the aerodynamic and IR characteristics of the nozzle require high-performance computing resources and advanced numerical methods.Simultaneously altering the values of all design factors of interest, while considering both aerodynamic and infrared characteristics, would considerably increase the computational workload.Therefore, most previous research on the elliptical double serpentine nozzle has focused on analyzing the influence of a single design factor or objective, which makes the conclusion limited.However, with the growing need for advanced stealth technology, there is an increasing demand for multi-objective design approaches that can simultaneously improve aerodynamic performance and reduce IR signatures [19,20].Conventional multi-objective optimization methods, such as Genetic Algorithms (GAs) require a considerably large number of objective function evaluations, rendering them impractical for this optimization problem due to the associated computational intensity [21,22].
This study aims to present a sample-efficient hybrid design optimization framework based on multi-objective Bayesian optimization (BO) to improve the infrared stealth characteristics of the serpentine nozzle while still emphasizing the enhancement of aerodynamic performance, which has significant engineering application potential and exploration value.The lengths of the first and second "S" sections and the aspect ratio of the nozzle exit are selected as the design variables.These are changed simultaneously during the iteration of the optimization.The rest of the paper is organized as follows: the parametric design of an elliptical double serpentine nozzle is first outlined in Section 2. A comprehensive optimization framework based on Bayesian optimization is presented in Section 3. A Bayesian-based multi-objective optimization framework is applied to optimize the shape of the elliptical double serpentine nozzle and the optimized models are compared with the baseline model in Section 4. Finally, the conclusions are summarized in Section 5.

Nozzle Design Methodology
The exhaust system under investigation comprises three primary sections: the exhaust mixer, the serpentine nozzle and the straight section.Based on a two-spool mixed flow turbofan engine at the design point, the exhaust mixer with two ducts is established.The hot gases flow through the core, while the low-temperature fluids flow through the bypass.Notably, the elliptical serpentine nozzle improves the stealth performance of the exhaust system in two ways.Firstly, the internal components of the exhaust system are obscured by the nozzle wall, and therefore, the infrared radiation and radar echo received by the detectors are reduced.Secondly, the elliptical-shaped outlet enhances the mixing of external and internal airflow, reducing the high-temperature exhaust gas temperature and thus lowering its infrared radiation.Figure 1 provides a brief description of the geometric parameters of the elliptical exit double serpentine nozzle.

Core
Center Cone a e b e Figure 1.Schematic diagram of the exhaust system under examination.The exhaust system comprises three main components: the mixing chamber, the double serpentine nozzle, and the straight section.For detailed geometrical parameters, refer to Table 1.The double serpentine nozzle with an elliptical exit can be divided into two individual single serpentine nozzles by the key section at the first bend.The axial length of the first serpentine section is denoted as L 1 with a vertical offset indicated by ∆Y 1 .The axial length of the second serpentine nozzle is L 2 with vertical offset ∆Y 2 .Then, the offset of the double serpentine nozzle can be expressed as ∆Y 2 − ∆Y 1 .Interfacing with the exit of the second serpentine nozzle is the straight section, which can restrict the gas flow and prevent nonuniform flow at the exit, with the length being L 3 .Detailed geometric parameters of the elliptical double serpentine nozzle are shown in Table 1.
Lee [23] proposed three distribution types which have superior aerodynamic performance in the subsonic diffusers.The so-called Lee curve is essentially a series of uniformly changing polynomial curves regulated by Equation (1), which is characterized by a slope of zero at both ends as shown in Figure 2. Given this characteristic, the distribution of area and curvature along the streamwise direction can also be governed by the Lee curve, as defined in Equations ( 2) and (3), respectively.The secondary flow inside the nozzle is influenced by the transverse pressure gradient, which is controlled by the centerline shapes.The streamwise pressure gradient that caused the flow separation is controlled by area distributions.The curvature distribution regulates the rate of transition from the circular inlet to the elliptical exit.In Equations ( 2)-( 4), L denotes the axial length of the nozzle and x refers to the location of the center along the axis.A i and A e represent the nozzle inlet area and exit area, respectively.κ e and κ i stand for the curvature of the nozzle exit and inlet, respectively.∆Y is the offset of the serpentine nozzle and y signifies the offset of the i-th section corresponding to the x coordinate.To reduce the loss instigated by flow separation at the first inflection section, the distribution of the centerline of the first serpentine section and the second serpentine section rapidly changes at the entrance and rapidly changes at the exit, respectively.In addition, a moderate area distribution is applied to acquire the section area and curvature for the serpentine nozzle.The specific coefficients for the Lee curve for the area, curvature, and centerline distribution are presented in Table 2.

First Section Second Section
Area,

Parametrization of Elliptical Serpentine Nozzle
A curvature-based control methodology is employed to generate the sections transitioning from a circular inlet to an elliptical exit.Given a specific nozzle inlet area A in , nozzle exit area A out , and aspect ratio AR of exit ellipse, the inlet and exit can be characterized by N discrete points, respectively.
In order to ensure the smoothness of the solid surface of the elliptical double serpentine nozzle generated in the 3-d construction, the discrete points of the curve imported into the CAD software must be evenly distributed.The uniform distribution of points on the elliptical nozzle exit is achieved by identifying the root of equal arcs length equation for the angle Φ.
To facilitate the transition from circular inlet to the elliptical exit, the curvature κ of inlet and outlet must first be resolved as follows: where s arc corresponds to the arc length divided by N points.In this study, N is set to 2000.Subsequently, the curvature of the j-th cross section is represented as: where κ in and κ out denote the curvature of the nozzle inlet and exit, respectively, and c j refers to the curvature distribution controlled by Equation (1).In this study, Equation ( 4) with coefficients of A = 3, B = −2, C = 0 and A = 4, B = −4, C = 1 is utilized to govern the curvature distribution of the first section and second section, respectively.It is worth noting that, due to truncation errors, the sum of the radians of each discrete arc length ∑ κ j • s may not equal 2π.This discrepancy implies that the cross section curve is not closed, which can prevent the successful generation of the nozzle wall.This issue can be addressed by applying a minor correction coefficient to each cross section.Once the curvature at each point along the cross section has been established, the coordinate of discrete points along the j-th curve can be iteratively calculated as follows: Following this, the coordinates of each cross-section are isotropically scaled to the desired area while retaining the same curvature.The area A for each cross-section is defined by Equation (4) with A = 3, B = −2, C = 0 for the two sections.Finally, based on the nozzle length L, nozzle offset ∆Y, and geometric parameters of the centerline, the final parametric characteristic curve of the serpentine nozzle can be derived by translating and rotating the section curves.

CFD Method
The outer mold lines around the nozzle are designed to represent an embedded propulsion system, and are aligned with the inflow boundary of the computational domain.A high-precision hexahedral structure grid is used to discretize the computational domain of the elliptical serpentine nozzle.The computational domain includes a rectangular far field bounding box containing free stream flow, with dimensions defined as 10D × 10D × 30D, where D is the diameter of the nozzle entrance, as shown in Figure 3a.Detailed grid configurations on the nozzle inlet and wall are shown in Figure 3b.Spalart-Allmaras (SA), standard k − ϵ and shear-stress transport (SST) k − ω turbulence model are used to carry out the numerical simulation for the elliptical double serpentine nozzle.Comparative analyses suggest that the SST k − ω model can accurately capture the characteristics of the pressure gradient flow and transonic shock wave with high convergence stability.In addition, many previous studies [1,15] have applied a SST k − ω turbulence model to the simulation of serpentine nozzles at different boundary conditions and obtained reasonable results.Therefore, the SST k − ω turbulence model in the Ansys FLUENT is adopted to simulate the flow characteristics of the elliptical double serpentine nozzle.The energy, momentum, density, turbulent kinetic energy, and turbulent dissipation rate are discretized using a second-order upwind scheme and the convergence accuracy of each equation is set to 1 × 10 −4 .The boundary conditions of the domains in the simulation are defined as follows: no-slip and adiabatic wall boundary conditions are imposed on the solid walls of the nozzle.The discrete ordinates (DO) radiation model is adopted to capture the radiation energies emitted from the solid wall.The external emissivity is set to 0.85 for the nozzle wall, cone and mixer.Pressure-inlet with total pressure and total temperature is applied to the core and bypass of the nozzle.The interior sections are meshed and set as the interior boundary, which is used to monitor the flow during the simulation.All simulations employ a nozzle pressure ratio (NPR) of 3.5 and an entrance temperature ratio of 0.4.The ambient domain is defined as pressure-far-field and the Reynolds number of the nozzle is 4.7 × 10 6 .The simulation is conducted at a flight altitude of 6 km and a cruising Mach number of 0.7, with the nozzle oriented at a zero angle of attack.The detailed values of the boundary conditions are presented in Table 3.A species transport model is selected to determine the spatial distribution of the absorption-emitting gas.The gas from the bypass section is considered to be air, whereas the exhaust gas from the core section is considered to be burned gas, with species properties detailed in Table 4.The arrangement of IR sensors is depicted in Figure 4. Sensors are strategically positioned on the xOy and yOz planes in a circular pattern, spaced at 2-degree intervals around the circumference and located 5 km from the center of the nozzle exit.In order to satisfy the requirements of the SST k − ω turbulence model, a boundary layer has been constructed near the solid wall inside the serpentine nozzle to ensure that y + < 5.As CFD simulations are time-intensive within the optimization process, solution convergence is accelerated by using the full multi-grid (FMG) initialization, which provides superior initial solutions.During the FMG initialization phase, the inviscid Euler equations are solved using first-order discretization to obtain an approximate solution.
To evaluate the grid size dependence of the numerical results, a grid-sensitive analysis for the serpentine nozzle is carried out.Four grids with different numbers of cells, where the total number of grids escalates from 2 million to 8 million, are simulated to characterize the aerodynamic characteristics of the serpentine nozzle.The ratio of static pressure at the exit to the total pressure at the nozzle inlet, denoted as p e /p * 0 , and the Mach number at the exit, represented as M a e , for different grids are extracted and exhibited in Figure 5.It is observed that the values of M a e and p e /p * 0 are almost indistinguishable between the medium mesh, dense mesh, and extremely dense mesh.Given that the infrared radiation calculation for the nozzle entails a highly complex ray tracing process, coupled with timeconsuming database queries and transmissivity calculations, an increase in the number of cells will significantly escalate the computation time.Therefore, to strike a balance between computational cost and precision, the mesh comprising 4.4 million cells has been selected to accurately capture the flow characteristics and temperature distribution of the elliptical double serpentine nozzle in the optimization process.It is worth noting that the CFD simulation for the grid with 4.4 million cells takes 8 h, whereas the infrared radiation calculation lasted 31 h.Consequently, the entire cycle for one nozzle model requires a total of 39 h.This computational task is executed on a 5.1 GHz Intel Raptor Lake i5-13600K CPU, equipped with 64 GB DDR5 Ram.Grid sensitivity analysis across four distinct grid sizes.To balance computational load with precision, the grid comprising 4.4 million cells was selected.

IR Radiation Calculation
The Reverse Monte Carlo ray tracing (RMCRT) method is adopted to calculate the infrared radiation of the serpentine nozzle.RMCRT is a stochastic method developed from the Monte Carlo Ray Tracing (MCRT) method, which can explicitly obtain solutions for the domain of interest without the need for the solution on the entire domain [24].Therefore, RMCRT is much faster than MCRT when simulating the IR signature of a large scene with a complex shape, like the exhaust system of an aircraft.Based on the three-dimensional flow field data obtained via CFD simulation, inclusive of the temperature, pressure and component concentration of absorbing-emitting gas on the face and voxel, RMCRT can calculate the IR radiation characteristics of the target via ray tracing.
Generally, the IR sensor is positioned at a considerable distance from the aircraft, which significantly exceeds the aircraft's characteristic length.This allows the aircraft to be treated as a point object emitting radiation.The radiant intensity emitted from the aircraft can be calculated as follows: where R is the distance from sensor to the target and λ 1 and λ 2 denote the boundaries of the detection waveband.E λ signifies the spectral irradiance directed toward the sensor, which can be determined by Unlike the traditional and widely adopted discrete transfer method (DTM) [25,26], the RMCRT method does not discretize the incident solid angle.Instead, it randomly emits N rays within the stereo angle space, employs a probability density function to determine the final absorption location of the rays, and substitutes the actual radiation target with N statistically significant radiation sources.Consequently, the irradiance can be written as where L bλ (i) denotes the black body spectral radiance at the absorption point of the i-th random ray, where θ i is the angle between the ray and the normal direction of the detection surface.Since the sensor is far away from the target, cos θ i ≈ 1, Equation (10) can be written as In order to solve Equation (11) and calculate the radiance over a certain spectral band of wavelength, the most precise approach to this involves looping through the waveband of interest and calculating the radiance at each wave point.This requires calculating the contribution of each spectral line for all molecules in the atmospheric layer; this is called a line-by-line (LBL) calculation.Although it is acknowledged as the most accurate model for simulating thermal radiative transfer, the LBL calculation is exceptionally time-consuming due to the small spectral resolution [27].To reduce the computation time, the statistic narrow band (SNB) model [28][29][30] is developed, wherein the spectral resolution is generally set between 5 and 50 cm −1 .Subsequently, Equation ( 11) can be simplified as where ∆λ is the spectral bandwidth and M denotes the number of discrete points in the waveband.
For the purpose of verifying the precision of the IR signature calculation code, the IR radiation characteristics of an experimental nozzle model proposed by Luo [31], as depicted in Figure 6 have been calculated and compared.The experimental nozzle is a scaled circular nozzle with a rear cone, making it suitable for comparison and validation.The simulation is conducted with parameters that mirror the experimental setup, where IR sensors are positioned 15 m from the nozzle exit in a circular arrangement.A comparison of experimental data and outcomes derived from the RMCRT is depicted in Figure 6.The experimental and simulated results only have a minor deviation at a few observation angles.For instance, at an observation angle of 20 • , the maximum deviation of 7.46% can be observed.Despite these minor discrepancies, the overall agreement between the results acquired from the reversed Monte Carlo ray tracing method and the experimental data demonstrates the reliability of the RMCRT method in predicting the IR radiation signature of the elliptical double serpentine nozzle.

Optimization Method and Process
The comprehensive optimization of the aerodynamic characteristic and infrared radiation signature of the elliptical double serpentine nozzle is a multi-disciplinary and multi-objective optimization problem.In this section, the definition of Pareto concepts involved in multi-objective optimization is restated.Then, the general framework of Bayesian Optimization (BO) is presented.Moreover, the key driver of adaptive sampling in multiobjective Bayesian optimization, acquisition function, is introduced.Finally, the description of the optimization problem of this study and a detailed flowchart of the optimization process are presented.

Background on Multi-Objective Optimization
In general, there are conflicts between objectives in most multi-objective optimization tasks.In the Pareto method, we say a solution is Pareto-optimal if it is not possible to move from that point and improve at least one objective function without having a detrimental effect on any other objective function.Herein, we consider minimizing M ≥ 2 objective functions f m (x) : X → R for m = 1, 2, . . ., M where X ∈ R xdim is the design space domain and xdim is the number of design parameters.We say a solution Then, the Pareto frontier can be defined as and the set of Pareto-optimal designs can be denoted as

General Framework of Bayesian Optimization
During the past three decades, a large number of multi-objective optimization algorithms have been proposed, such as Evolutionary Algorithms (EA, e.g., NSGA-II) and Swarm Intelligence (SI, e.g., MOPSO), which show strong robustness to many problems.However, these algorithms require a large number of function evaluations [21,22], rendering them infeasible for solving the optimal problems with expensive-to-evaluate "black-box" functions (e.g., CFD simulation, hyper-parameter tuning in Machine Learning, etc.).The shape optimization of an elliptical double serpentine nozzle in this study involves numerous CFD simulations and infrared radiation calculations, which is suitable to be solved under multi-objective Bayesian optimization (MOBO).MOBO offers a sample-efficient method by generating the candidate solutions that best balance the predicted value and the uncertainty given all previous data [32].
The general framework of MOBO is shown in Figure 7.A probabilistic surrogate model is initially constructed once the initial sampling is completed by acquiring the corresponding objective values by implementing the black-box function.Typically, the surrogate is a Gaussian process regression (GPR), which is a flexible, non-parametric model known for its well-calibrated predictive uncertainty [33].Subsequently, based on the posterior distribution P ( f |D) over the true function values f provided by the GPR, the probabilistic surrogate model is updated via adaptive sampling, where the infill design point is determined by finding the maximum of the acquisition function x * = arg max x∈X α(x).This process is repeated for a fixed number of iterations, determined by the computational budget.Once the BO is finished, the probabilistic model is optimized to obtain the Pareto set of the optimization problem.The key to the success of BO is that the acquisition function is cheaper to evaluate than the black-box function [34].In this study, the black-box function is the CFD simulation and IR calculation of the nozzle.Therefore, this can save a lot of expensive evaluations.

Expected Hyper-Volume Improvement
During each iteration of MOBO, an infill point is selected to enhance the surrogate model's accuracy.This selection is based on optimizing a predifined acquisition function.The expected improvement (EI) criterion is a widely used acquisition function in Bayesian optimization, quantifying the potential improvement of the surrogate model at a specific input location.There are several variants of acquisition functions based on EI, such as Euclidean-based EI (EEI) [35], expected hypervolume improvement (EHVI) [36], ParEGO [37], and expected inverted penalty boundary intersection improvement (EIP-BII) [38].Among these variants, EHVI has been shown to perform well when applied to both artificial and aerodynamic problems, which makes it a highly promising MOBO method for multi-objective aerodynamic design optimization [39].
EHVI was first proposed by Emmerich [36] at 2006, and it is grounded in two fundamental concepts: hypervolume (HV) and hypervolume improvement (HVI).Both are utilized as general performance indicators in multi-objective optimization [40].The hypervolume indicator of a finite approximate Pareto frontier P is the M-dimensional Lebesgue measure λ M of the space dominated by P and bounded from below by a reference point r: where [r, υ] denotes the hyper-rectangle bounded by vertices r and υ.Typically, the ref- erence point r is specified by the decision maker [41].In this paper, the reference point was set to [1,1] after the objective functions were normalized using the minimum and maximum values from the current experimental design.The HVI of a set of points P ′ with reference to an existing approximate Pareto frontier P and reference point r is defined in Equation ( 16) and an illustration that depicts computing the HVI with new candidate f (x) is shown in Figure 8.
HVI(P ′ |P, r) = HV(P ∪ P ′ |r) − HV(P |r) EHVI is the expected value of how much hypervolume improvement (HVI) can be obtained considering the current non-dominated solutions, the reference point, the Kriging prediction, and the Kriging estimation error.The analytical expression of EHVI can be formulated as where M is the number of objectives and ϕ( f ) is the probability density function.
To verify the feasibility of the qEHVI acquisition function, we evaluate one synthetic benchmark, the Branin-Currin problem [42] (d = 2, M = 2), and one real-world benchmark, the Vehicle Safety problem [43] (d = 5, M = 3).The log hypervolume difference, which is defined as the difference between the hypervolume of the true Pareto front and the hypervolume of the approximate Pareto front based on the current observed data, is used as the evaluation criterion.The maximum hypervolume, reference point, search space and detailed problem definitions are presented in Appendix A. Four algorithms, including random search (Sobol), NSGA-II, and Bayesian optimization with qParEGO and qEHVI acquisition functions, respectively, are carried out on the two benchmark functions.In each run of the experiments, all algorithms are provided with 10 initial evaluations, and these initial inputs are shared across all algorithms.The experiments employed the Bayesian optimization model and acquisition function from the BoTorch framework [44], while the default implementations of the NSGA-II algorithm provided by the pymoo library [45].
Figure 9 illustrates the mean logHV difference over 20 runs.It can be seen that qEHVI outperforms the qParEGO, NSGA-II, and Sobol in terms of sequential optimization performance for the two optimization problems.Following the random initialization phase, the logHV difference of MOBO with qEHVI rapidly decreases.This indicates an increase in the HV and a closer proximity between the obtained Pareto frontier and the true frontier.The final experimental HV after 100 iterations and ideal HV for the two problems are shown in Table 5.It can be observed that qEHVI achieves a final HV that is very close to the theoretical HV.Compared to the Sobol, NSGA-II and qParEGO methods, qEHVI exhibits HV improvements of 38.84%, 2.69%, and 0.07% for the Branin-Currin problem, and 30.54%, 4.37%, and 0.264% for the vehicle safety problem, respectively.In addition, it is worth highlighting that the Bayesian algorithm demonstrates the lowest variance among the algorithms across multiple runs, which is particularly significant in practical applications where the expensive nature of the setting precludes the ability to perform multiple runs.This characteristic of the Bayesian algorithm makes it an attractive choice for the aerodynamic and stealth design, where computational resources are limited.
HV log_diff = log(HV max − HV cur ) .The illustration of computing the hypervolume improvement with infill point f (x) when there are two objectives.The current Pareto frontier P is given by the green points and the red area is the hypervolume surrounded by the P and given reference point r.The white rectangles S 1 , . . ., S K represent a disjoint box decomposition of the non-dominated space, which can be utilized to compute the hypervolume improvement efficiently.

Formalizing the Problem
The aerodynamic performance of the elliptical double serpentine nozzle can be evaluated using the thrust coefficient C f , which is defined as the ratio of the actual thrust of the nozzle F x to the ideal thrust F i of the one-dimensional isentropic fully expanded flow, as shown in Equation (19).The closer the thrust coefficient is to 1, the smaller the flow loss inside the nozzle and the better aerodynamic performance of the nozzle.
and the F x and F i can be formulated as where subscripts in and e represent the inlet and exit sections of the nozzle, respectively.ṁi is the mass flow rate, which can be formulated as In this study, the thrust coefficient of the serpentine nozzle C f is used as the objective function to evaluate the aerodynamic performance of the nozzle.
The objective function of the infrared radiation signature of the nozzle is a weighted summation of the infrared intensity mean values at different detection angles in different wavebands, as defined in Equation (23).The smaller the infrared intensity, the better the stealth performance of the elliptical double serpentine nozzle.
where N is the number of sensors, I is the infrared radiation intensity in the MWIR and LWIR band of the nozzle simulated by RMCRT.The primary objective of the elliptical double serpentine nozzle is to mitigate the IR signature by providing shielding for the high-temperature section and facilitating the thorough mixing of the high-temperature exhaust with the cool ambient air.However, this inevitably affects the aerodynamic performance of the engine exhaust system.Therefore, in order to search for the good infrared radiation suppression effect while taking into account the aerodynamic performance of the elliptical double serpentine nozzle, the following multi-objective optimization problem can be defined: where C f and f ir represent the thrust coefficient and infrared radiation objective function as defined in Equations ( 19) and ( 23), respectively.AR is the aspect ratio of the elliptical exit of the nozzle; L 1 is the length of the first serpentine section.

Optimization Procedure
Despite the existence of several space-filling techniques that are well suited for the design of experiment (DOE), Latin Hypercube Sampling (LHS) has attracted the most attention in recent years [46,47].The Latin hypercube sampling (LHS) technique allows designers to freely select the number of design parameters and observations, rendering it an appealing option for optimization problems involving computationally expensive objective functions, such as aerodynamic simulations.In this study, an optimal Latin hypercube sampling (OLHS) is adopted to generate initial samples within the search space.This approach enhances LHS's space-filling properties by utilizing a columnwise-pairwise algorithm to optimize the distribution of samples based on the S-optimality criterion [48,49].
The overall optimization procedure is shown in Figure 10.Firstly, 40 DOE points are generated as the initial samples in the search space.Subsequently, the aerodynamic characteristics and IR properties of 40 nozzles are then evaluated using the CFD and RMCRT method, respectively.Based on these data, a probabilistic surrogate model is constructed using the Kriging model to sequentially update the solution towards the Pareto frontier.At each iteration of MOBO, a new sample is generated according to the optimization of the acquisition function.This iterative process continues until the termination criteria are met.Finally, the non-dominated sets are obtained, from which the optimal serpentine nozzles are selected.

Optimization Details
Following the OLHS DOE process, a total of 40 samples are obtained, and both the aerodynamic characteristics and IR signature are simulated for each sample.The optimization objective functions, as defined in Equations ( 19) and ( 23), are subsequently computed, forming the initial dataset with which the probabilistic surrogate model is constructed.Prior to the optimization process, the Kriging surrogate model is assessed for its ability to generalize given data.The leave-one-out cross-validation (LOOCV) method is adopted to evaluate the Kriging surrogate model's ability to predict objective values at unobserved locations.A detailed description of the LOOCV procedure can be found in [50].
The diagnostic plots from the LOOCV procedure for the C f and f ir surrogate are presented on the left side of Figure where values on the horizontal axis denote the exact results obtained through numerical simulation, and values on the vertical axis represent the prediction results from the Kriging model when the corresponding sample is omitted.A high-performing surrogate model would be characterized by samples located in close proximity to the diagonal line.In other words, the closer the samples are to the diagonal line, the better the surrogate model's performance when it comes to capturing the underlying patterns and accurately predicting objective values at unobserved locations.It can be seen that the majority of observations cluster closely around the diagonal line, with only a few observations slightly deviating from the ideal location, suggesting that the model is capable of providing reliable predictions for both the thrust coefficient and IR signature of the elliptical double serpentine nozzle.The right panel of Figure 11 depicts the distributions of the cross-validated residuals.The standardized residuals within a three-sigma range, denoted by the red shaded areas, correspond to a 99.7% confidence interval for Gaussian processes like Kriging [51].For the aerodynamic Kriging model, all samples fall within the desired interval, signifying the robust performance of the constructed surrogate model in predicting the thrust coefficient of the elliptical double serpentine nozzle.With respect to the IR surrogate model, three observations marginally exceed the desired boundaries.Nevertheless, apart from these anomalies, the majority of the residuals of the observations are distributed around the zero standard, indicating the model's sufficiency in accurately predicting the IR signature of the serpentine nozzle.
In summation, the findings confirm the sufficient accuracy of the surrogate models constructed to predict both the aerodynamic performance and the IR characteristics of the elliptical double serpentine nozzle.These models, therefore, represent valuable tools for optimizing the nozzle's design parameters and assessing their impact on the thrust coefficient and IR signature.This, in turn, contributes to the advancement of high-performance nozzle development for a range of applications.
In this study, due to the high computational cost of the CFD and IR simulations, only one new candidate is generated per iteration (q = 1) of Bayesian optimization.A straightforward initialization heuristic is employed to select 20 restart initial locations from a pool of 256 randomly generated points.A multi-start optimization of the qEHVI acquisition function is undertaken using LBFGS-B, with exact gradients computed via auto-differentiation powered by PyTorch.
Figure presents the entire set of solutions projected into the infrared intensity and C f space.This includes both the initial samples and candidate observations accumulated throughout the iterations of Bayesian optimization.Since the multi-objective optimization problem defined above is a minimization problem, the x-axis denotes the negative C f .The brightness of the solutions' color increases as the iterations progress in the BO.Such visualization of the results, preferably including interactive features, is particularly important to analyze in multi-objective optimization, as it can provide valuable insights for designers or decision-makers when evaluating different design option.It can be observed that the whole set of candidates obtained using the Bayesian algorithm appears to cluster in a more favorable area (the lower left corner of the Figure 12), which reflects the advantage of incorporating prior domain knowledge into Bayesian algorithms.This finding, together with the previous observations made on the respective performance of the algorithms in Figure 9, suggests that the Bayesian algorithm is capable of finding higher-quality trade-offs among the aerodynamic performance and infrared stealth performance of the elliptical double serpentine nozzle.

Comparison of Baseline and Optimal Nozzles
Furthermore, in order to highlight the characteristics of the elliptical double serpentine nozzle, the simulation results of an ordinary single serpentine nozzle with an elliptical exit are also presented for comparison.The geometric details of these models are shown in Figure 13.All four models have the same vertical offset.Compared to the Optimal-1 and Optimal-2 model, the first bend length of the Baseline model is smaller.Notably, the Optimal-2 model has a significantly larger aspect ratio than the others.
To highlight the superiority of each model, a thorough analysis of the aerodynamic and infrared radiation characteristics for each model, as presented in Figure 13, is conducted.The comparative simulation outcomes, detailed in Table 6, reveal that the Optimal-1 model exhibits enhanced aerodynamic characteristics and infrared stealth performance compared to the Baseline model.Specifically, the thrust coefficient C f improves by 0.27% and the infrared radiation intensity f IR declines by 10.08%.In contrast, the thrust coefficient of the Optimal-2 model has deteriorated due to its considerable aspect ratio.Nevertheless, this model achieves a more significant reduction in IR intensity than Optimal-1, with an average infrared intensity reduction of 18.12% across MWIR and LWIR.This reduction is attributed to the larger aspect ratio of the nozzle exit, which better shields the engine's hot components and enhances the blending of the hot plume with ambient air.The aerodynamic performance of the single serpentine nozzle improves by 0.61% compared with the Baseline model.However, since the nozzle wall cannot fully block the hot part inside the nozzle, the mean IR intensity has been increased by 9.36%.Furthermore, a substantial IR intensity peak is observed at ϕ = 80 • , where the hot components such as the turbine section are directly exposed to the IR sensor, as depicted in Figure 14.It is noteworthy that while the double serpentine nozzle can detrimentally impact the aerodynamic performance of the nozzle flow field, the complete shielding of the high-temperature part and hot plume inside the nozzle can effectively reduce the infrared radiation signature.The data presented in Table 6 suggest that there is a trade-off between aerodynamic performance and infrared stealth performance for the elliptical double serpentine nozzle.By optimizing the nozzle shape, it is feasible to achieve an overall performance enhancement relative to the Baseline model.
The z-axis time-averaged vorticity distribution across different cross sections of each model is presented in Figure 15, allowing us to analyze the internal flow field within the nozzle.As illustrated, the presence of the central cone results in an augmented rotational velocity of the time-averaged vorticity as the radius of the cone decreases, in accordance with the law of conservation of angular momentum.Consequently, a concentrated vortex is formed downstream of the central cone.The vortex generated by the central cone then accelerates the expansion of the core flow, thereby enhancing the mixing of the core and bypass flows.Furthermore, the shear effect between the core and bypass flows induces a pronounced streamwise vortex within the elliptical nozzle.The intensity of these streamwise vortexes is notably higher in the double serpentine nozzle compared to the single serpentine nozzle due to its curved geometric configuration.In the case of the Baseline model, where the length of the first bend is shorter than the other two double serpentine models, a substantial flow deflection occurs within this section.Consequently, two distinct streamwise vortexes develop within the first bend and persist until they reach the exit of the nozzle.In contrast, the Optimal-1 and Optimal-2 models, which have longer first bend sections, develop two pairs of streamwise vortexes in the second bend of the nozzle.As a result, the flow at the nozzle exit becomes more uniform, and the aerodynamic performance of the nozzle improves.The IR radiation intensity distributions of the single serpentine nozzle and Baseline model are illustrated in Figure 14.The polar lines on the left represent the total IR radiation intensity, which is the summation of the IR radiation of the solid wall of the nozzle (right) and the IR radiation of gas (middle).The IR radiation intensity of the wall received by the detector is the radiation emitted from the solid wall boundaries of the nozzle, including the nozzle wall, the nozzle mixer, and the nozzle central cone, which has been attenuated by the gas.The IR radiation intensity from the nozzle wall is determined by the visibility of the nozzle's hot part and the wall itself to the detector.Given that the IR radiative power E is significantly sensitive to the temperature (E = σT 4 ), the temperature distribution of the visible wall is the primary determinant of the nozzle's IR radiation characteristic.This sensitivity results in a noticeable peak in the IR intensity for the single serpentine nozzle, as the single serpentine nozzle does not fully block the hot part inside the nozzle.
The primary effect of the double serpentine nozzle is to shield the high-temperature components inside the nozzle, such as the nozzle central cone and turbo section.The IR intensity of the Baseline model remains relatively low at most detection angles, with the exception of when the detection angle θ falls within the range of 90 • -140 • .This exception is attributable to the hot streak located at the upper nozzle wall, a phenomenon which will be further explained in the following discussions.The hot exhausted gas produced by the engine can be divided into two parts: gas inside the nozzle exit and gas outside the nozzle exit.The single serpentine configuration would require a significant offset for complete internal component shielding, a requirement that presents substantial design and integration challenges, as noted by Cheng et al. [1].
The top panel of Figure 16 displays the temperature contours across the cross sections of the Baseline model, Optimal-1 and Optimal-2 models.As the flow develops, the hightemperature core flow consistently mixes with the bypass flow, gradually reducing the temperature of the core flow.The wall curvature and the transition from circular to elliptical induce the formation of two vortices at the sides of the nozzle.This leads to a distinct temperature distribution pattern at the elliptical exit of the elliptical double serpentine nozzle.Notably, the non-uniformity at the nozzle exit of Optimal model is appreciably diminished.
The temperature distributions of each model's nozzle wall are depicted in the bottom panel of Figure 16.Notably, a band-like distribution of localized high-temperature regions, also known as a "hot streak", is evident on the upper wall of the nozzle.In the Baseline model, this hot streak originates at the second bend of the nozzle and extends into the straight sections.Conversely, the hot streak in the Optimal models is primarily confined to the straight sections.In comparison to the Baseline model, not only is the area of these hot streaks smaller, but their temperature is also significantly reduced.The temperature manifests a gradual decrease from the central region towards the side in the cross-sectional direction of the nozzle.Figure 17 depicts the angular distribution of the infrared intensity of the elliptical double serpentine nozzle on the pitch detection plane.The top panel represents the infrared radiation in the middle-wavelength band (3-5 µm), while the bottom panel corresponds to the infrared radiation in the long-wavelength band (8-14 µm).There is a clear distinction between the distribution patterns of the MWIR intensity and those in the LWIR.In the MWIR range, the infrared radiation is dominated by the hot plume of the nozzle, while the solid radiation accounts for a relatively small part.As the detector shifts away from the aft end of the nozzle (θ = 90 • ), the MWIR intensity gradually increases, as demonstrated in Figure 17b.This rise is attributable to the increasing projected area of the hot plume on the detector focal plane as the detector moves away, thereby increasing the infrared intensity.When the detector is positioned directly above (θ = 0 • ) or below (θ = 180 • ) the nozzle, only the infrared radiation emitted from the plume is accessible to the detector.
It can be observed that the hot streak has a significant impact on the infrared radiation characteristic of the three elliptical double serpentine nozzles.The peak solid infrared radiation within both the MWIR (Figure 17c) and LWIR (Figure 17f) occurs at a detection angle around 110 • in the vertical detection plane.This is because the infrared radiation emitted from the hot streak can be received by the IR sensors, and the projected area is maximized at this detection angle.In addition, the solid LWIR intensity is larger than the MWIR intensity, which can be explained by Wein's displacement law.Given the area and temperature of the hot streak in Optimal-1 and Optimal-2 model, these models exhibit a noticeable decrease in solid infrared radiation, which is consistent with the temperature distribution pattern of the hot streak shown in Figure 16.
Observing Figure 17b, it is noticeable that in the middle waveband, the gas radiation of the Optimal-2 model (blue line) is smaller than the gas radiation of the Optimal-1 model (gray line) at most detection angles.This reduction is attributed to the larger exit aspect ratio of the Optimal-2 model, which enhances the mixing of the hot exhaust flame with ambient cool air, thereby reducing the temperature of the tail flame core area (as shown in Figure 18) and consequently decreasing its radiation energy.On the other hand, the larger aspect ratio of Optimal-2 model results in a somewhat stronger hot spot in the Optimal-2 model compared to the Optimal-1 model.Figure 17f reveals that the solid radiation of the Optimal-2 model is slightly higher than that of the Optimal-1 model.Simultaneously, the stronger vortex caused by the larger aspect ratio in the Optimal-2 nozzle leads to a minor reduction in the aerodynamic characteristics compared to the Baseline model and the Optimal-1 model.Moreover, in our study, we found that for the serpentine nozzle, due to the complete shielding of high-temperature components by the double S-bend design of the nozzle, the solid radiation emanates predominantly from the nozzle wall.At most detection angles, the radiation from other components is virtually negligible.This finding underscores the effectiveness of the double S-bend nozzle in protecting against infrared radiation from the internal high-temperature sections.Figure 18 reveals the infrared radiation intensity images within the MWIR and LWIR bands for different models.These images provide an intuitive view of the variations in high-temperature areas of different nozzles.The left side of Figure 18 presents the infrared radiation intensity images observed from above the nozzle (θ = 90 • , ϕ = 0 • ).It is evident that there is an alternation between high-and low-temperature zones at the nozzle exit generated by shock waves.The Optimal-1 model exhibits a high-temperature tail flame length nearly identical to that of the Baseline model.However, the exhaust flame length of the Optimal-2 model is noticeably shorter, with bifurcation evident in the plume.This is because the Optimal-2 model's nozzle exit has a larger aspect ratio (AR = 7.8), which enhances the mixing of the plume with the ambient cold air, resulting in the strong 3-D effect of the fluid.Thereby, the high-temperature core length of Optimal-2 is shorter than that of the other two models.Correspondingly, the IR radiation intensity of gas is also reduced, as shown in Figure 17.The right side of Figure 18 features the infrared radiation intensity images from the rear perspective (θ = 90 • , ϕ = 90 • ).The Baseline model exhibits a significantly higher hot streak intensity compared to the Optimal-1 and Optimal-2 models, resulting in a notably dimmer appearance for the latter.This observation corroborates the quantitative findings presented earlier in Figure 17 and Table 6.

Conclusions
In this paper, a state-of-the-art multi-objective Bayesian optimization algorithm was adapted to an aerodynamic and infrared stealth optimization problem for the elliptical double serpentine nozzle.Based on the curvature control method, a series of nozzles were designed to perform the Bayesian optimization.The Pareto frontier was finally obtained and two non-dominated solutions on the PF were selected to obtain the optimal model.A detailed comparison of the characteristics of these two models with the Baseline model and one single serpentine nozzle was carried out.The main conclusions can be drawn as follows: • The length of the first bend and the aspect ratio of the elliptical double serpentine nozzle have an important influence on the performance of the nozzle.Notably, the Optimal-1 model demonstrated enhanced infrared radiation and aerodynamic performance following the optimization, with a 10.08% decrease in infrared radiation intensity and a 0.27% improvement in the thrust coefficient.

•
The Optimal-2 model, with a larger aspect ratio, showed an impressive improvement in the infrared radiation performance, demonstrated by an 18.12% reduction in infrared radiation intensity.However, the larger transition from circular inlet to elliptical exit also deteriorated the flow field of the Optimal-2 model, while the thrust coefficient decreased by 0.14%.

•
The optimization framework introduced in this paper effectively tackles the complex multi-objective optimization problems of the elliptical double serpentine nozzle.This framework not only provides reliable technical guidance for future research, but is also a promising approach to complex and costly optimization challenges in aircraft design, including the improvement of aerodynamic, radar, and infrared stealth capabilities. where

Figure 3 .
Figure 3. Computational grid of the elliptical double serpentine nozzle.(a) global computational grid.(b) nozzle grid.

Figure 4 .
Figure 4. Arrangement of infrared sensors.IR sensors are strategically positioned on the xOy and yOz planes in a circumferential pattern.They are spaced at 2-degree intervals around the circumference, with each sensor located 5 km from the nozzle.

Figure 5 .
Figure5.Grid sensitivity analysis across four distinct grid sizes.To balance computational load with precision, the grid comprising 4.4 million cells was selected.

Figure 6 .
Figure 6.Verification of RMCRT.Panel (a) Model of the experimental nozzle depicts a scaled model with a rear cone proposed by Luo [31].Panel (b) Comparison of the simulation result and experimental data for IR validation, validates the accuracy of the RMCRT method in simulating the infrared radiation signature of the nozzle.

Figure 7 .
Figure 7. General framework of multi-objective Bayesian optimization.

Figure 8
Figure 8.The illustration of computing the hypervolume improvement with infill point f (x) when there are two objectives.The current Pareto frontier P is given by the green points and the red area is the hypervolume surrounded by the P and given reference point r.The white rectangles S 1 , . . ., S K represent a disjoint box decomposition of the non-dominated space, which can be utilized to compute the hypervolume improvement efficiently.

Figure 9 .
Figure 9.A comparison of different algorithms used for synthetic test function (left) and a realworld problem (right).The hypervolume indicator is shown with regard to the number of function evaluations.Every experiment has a batch size of 4 and 25 iterations in total.The means and one standard error across 10 trials are plotted.

Figure 11 .
Figure 11.Validation of the Kriging surrogate models.Top: cross validation of C f surrogate model.Bottom: cross validation of f ir surrogate model.(Left) LOOCV plots for the C f and f ir Kriging models.The dashed diagonal line represents the locations of perfect predictions.(Right) Plots of standardized cross-validated residuals for C f and f ir Kriging models.Red region indicates the interval of 3 standard deviations.

Figure 12 .
Figure 12.Collection of MOBO solutions projected the objective space where the color corresponds to the MOBO iteration at which the point was collected.It includes an initial set of 40 sample points, with an additional 20 points derived via the optimization of the acquisition function.The red line delineates the Pareto frontier achieved through iterative optimization.The baseline model is represented by the triangular marker.Two optimal models are chosen for in-depth analysis: one excelling in aerodynamic performance and the other in infrared stealth performance.

Figure 13 .Figure 14 .
Figure 13.Geometric outlines of the baseline model, single serpentine nozzle, Optimal-1, and Optimal-2 model.The top panel presents the side view of the nozzles, the middle panel illustrates the exit profiles of the nozzles, and the bottom panel provides an isometric perspective of the nozzles.

Figure 17 .
Figure 17.Comparison of infrared intensity between each model.(a) Total infrared radiation intensity in MWIR range.(b) Gas infrared radiation intensity in MWIR range.(c) Solid infrared radiation intensity in MWIR range.(d) Total infrared radiation intensity in LWIR range.(e) Gas infrared radiation intensity in LWIR range.(f) Solid infrared radiation intensity in LWIR range.

Figure 18 .
Figure 18.Comparison of infrared radiance image between each model.(Left) Observed from above the nozzle.(Right) Observed from directly behind the nozzle.

Table 1 .
Design parameters of the elliptical serpentine nozzle.

Table 2 .
Distribution parameters of the elliptical serpentine nozzle.

Table 4 .
Properties of the exhaust gas (mass fraction).

Table 5 .
Mean hypervolume across 20 runs per algorithm at the final evaluation.

Table 6 .
Parameters and results of each model.