Moritz Krügener
email: moritz.kruegener@gmail.com

Jose Felix

Zapata Usandivaras
email: jose.zapata-usandivaras@isae-supaero.fr

Michaël Bauerheim
email: michael.bauerheim@isae-supaero.fr

Annafederica Urbano
email: annafederica.urbano@isae-supaero.fr

Coaxial-Injector Surrogate Modeling Based on Reynolds-Averaged Navier-Stokes Simulations Using Deep Learning

published or not. The documents may come

Coaxial-Injector Surrogate Modeling based on Reynolds-Averaged Navier-Stokes Simulations using Deep Learning

M. Krügener * , J. F. Zapata Usandivaras † , M. Bauerheim ‡ , A.

I. Introduction

Liquid Rocket Engines (LRE) are complex systems that are composed of several sub-systems (cooling system, feed system, thrust chamber...) that are designed to be reliable and optimised in order to obtain the highest performance and the lowest mass. The process is iterative, and in the preliminary sizing phase, heavily relies on low order semi empirical correlations calibrated on existing engines. These low fidelity models induce high uncertainties in the subsystems designs that propagate and affect the global engine sizing. As a consequence, the engine development strongly relies on experimental studies, that go up to hot-tests of the engine, making use of very expensive experimental facilities. Increasing the numerical modelling reliability is a strategic matter that needs to be addressed in order to reduce the development cost of innovative rocket engines for future launchers, in particular LRE, and thus ensuring access to space. The combustion chamber and injector design are key elements of any rocket engine [START_REF] Martin | Flame Characteristics of a High-Pressure LOX/H2 Rocket Combustor With Large Optical Access[END_REF]. High fidelity numerical simulations, and in particular Large Eddy Simulations (LES), have demonstrated their ability to finely describe complex physical phenomena in LRE combustion chambers [START_REF] Urbano | Exploration of Combustion Instability Triggering Using Large Eddy Simulation of a Multiple Injector Liquid Rocket Engine[END_REF][START_REF] Urbano | Driving and Damping Mechanisms for Transverse Combustion Instabilities in Liquid Rocket Engines[END_REF][START_REF] Maestro | Large Eddy Simulation of Combustion and Heat Transfer in a Single Element GCH4/GOx Rocket Combustor[END_REF][START_REF] Huo | Large-Eddy Simulation of Supercritical Combustion: Model Validation Against Gaseous H2-O2 Injector[END_REF]. Nevertheless, LES is still a high computational cost methodology, with long computational times, which can limit its use in the context of concrete applications. An alternative, when only stationary operating conditions are of interest, is to use Reynolds-Averaged Navier-Stokes (RANS) simulations [START_REF] Riedmann | Modeling of Combustion and Flow in a Single Element GH2/GO2 Combustor[END_REF][START_REF] Perakis | Heat Flux and Temperature Evaluation in a Rectangular Multi-Element GOX/GCH4 Combustion Chamber Using an Inverse Heat Conduction Method[END_REF][START_REF] Lapenna | Modeling the Equations of State Using a Flamelet Approach in LRE-Like Conditions[END_REF][START_REF] Perakis | Experimental and Numerical Investigation of CH4/O2 Rocket Combustors[END_REF][START_REF] Lempke | Unsteady High-Order Simulation of a Liquid Oxygen/Gaseous Hydrogen Rocket Combustor[END_REF]. However, even restitution times of RANS are too high if used in an optimisation framework loop of a whole engine. What is the best way to exploit the capability of Computational Fluid Dynamics (CFD) in this context is an interesting question which is under investigation today [START_REF] Wang | A High-Fidelity Design Methodology Using LES-Based Simulation and POD-Based Emulation: A Case Study of Swirl Injectors[END_REF][START_REF] Mak | An Efficient Surrogate Model for Emulation and Physics Extraction of Large Eddy Simulations[END_REF]. In order to answer this question, a strategy is proposed in the present paper using models of various fidelity. The main idea is to build surrogate models with real-time inference capabilities for the injector and the combustion chamber, using machine learning algorithms and data from numerical simulations. The obtained surrogate models can then be used in the framework of an optimisation problem [START_REF] Gray | Chemical-Equilibrium Analysis with Adjoint Derivatives for Propulsion Cycle Analysis[END_REF][START_REF] Chen | Numerical Investigation of Minimum Drag Profiles in Laminar Flow Using Deep Learning Surrogates[END_REF]. Once deployed, the surrogate model will allow a large number of different configurations to be spanned, providing the results of the CFD simulations with a very low inference time, of the order of the millisecond. In the context of design and optimisation of the whole rocket engine, the surrogate model of the chamber can permit to increase the accuracy of the final configuration performance. This eventually can reduce the number of experimental tests and consequently the whole development cost. The derivation and use of data driven surrogate models in the aerospace community is an active field of research. In particular, response surface models have been applied in a wide range of applications including nozzle shape optimization [START_REF] Simpson | Kriging Models for Global Approximation in Simulation-Based Multidisciplinary Design Optimization[END_REF] and rocket injector design [START_REF] Shyy | Response Surface and Neural Network Techniques for Rocket Engine Injector Optimization[END_REF][START_REF] Vaidyanathan | Computational-Fluid-Dynamics-Based Design Optimization for Single-Element Rocket Injector[END_REF]. In these works, surrogate models are either interpolators, usually high polynomials surface functions, or regressors, usually based on gaussian processes. A survey of these methods can be find for instance in the paper from Hwang and Martins [START_REF] Hwang | A Fast-Prediction Surrogate Model for Large Datasets[END_REF]. Recently, the great developments of artificial intelligence algorithms have been recognised in the fluid mechanics community and their application to CFD is a growing field of investigation [START_REF] Brunton | Machine Learning for Fluid Mechanics[END_REF]. In particular, deep learning (DL) has shown a great potential for both modelling and analysing of complex systems, including 2D and 3D geometries [START_REF] Baqué | Geodesic Convolutional Shape Optimization[END_REF][START_REF] Tompson | Accelerating Eulerian Fluid Simulation With Convolutional Networks[END_REF], flow control and optimisation [START_REF] Chen | Numerical Investigation of Minimum Drag Profiles in Laminar Flow Using Deep Learning Surrogates[END_REF]. The use of deep learning for regression on CFD data, in order to obtain data based models, has been investigated in the framework of both steady external aerodynamics [START_REF] Chen | Numerical Investigation of Minimum Drag Profiles in Laminar Flow Using Deep Learning Surrogates[END_REF], unsteady fluid mechanics [START_REF] Ajuria | Towards an Hybrid Computational Strategy Based on Deep Learning for Incompressible Flows[END_REF] , acoustics [START_REF] Alguacil | Predicting the Propagation of Acoustic Waves Using Deep Convolutional Neural Networks[END_REF], and internal reactive flows problems [START_REF] Ma | Supervised Learning Mixing Characteristics of Film Cooling in a Rocket Combustor Using Convolutional Neural Networks[END_REF][START_REF] Shadram | Neural Network Flame Closure for a Turbulent Combustor with Unsteady Pressure[END_REF]. Of particular interest is the work by Thuerey et al [START_REF] Thuerey | Deep Learning Methods for Reynolds-Averaged Navier-Stokes Simulations of Airfoil Flows[END_REF], who investigated the accuracy of deep learning for the inference of the flow around an airfoil. The developed surrogate models, based on U-net, a particular class of Convolutional Neural Networks, were able to reproduce the pressure and velocities fields with errors smaller than 3%. Despite the increasing number of works in this domain, the use of these approaches applied to rocket engines remains limited [START_REF] Ma | Supervised Learning Mixing Characteristics of Film Cooling in a Rocket Combustor Using Convolutional Neural Networks[END_REF][START_REF] Shadram | Neural Network Flame Closure for a Turbulent Combustor with Unsteady Pressure[END_REF][START_REF] Mcquarrie | Data-Driven Reduced-Order Models Via Regularised Operator Inference for a Single-Injector Combustion Process[END_REF][START_REF] Haemisch | LUMEN -Design of the Regenerative Cooling System for an Expander Bleed Cycle Engine Using Methane[END_REF] and to the authors knowledge there are no work available focusing on their application for injector design.

The aim of the present paper is to carry out a proof of concept of the proposed methodology based on machine learning by applying it on a simple configuration. The selected test case is a single coaxial injector chamber, using gaseous methane and oxygen as propellants. The database will be generated using 2D axisymmetric RANS simulations.

Specific geometric parameters (for instance the chamber diameter and injector diameters) and specific operative conditions are defined as design parameters. The objective is to build surrogate models able to provide, with a very low inference time, three different types of quantities: performance parameters (0D), the wall heat flux distribution (1D) and a visualisation of the chamber temperature field (2D). In order to achieve this objective several research questions will be addressed. First, the paper will investigate how to deploy a strategy in order to build a design of experiments from the selected case. The second question addresses the type and architecture of the machine learning tool necessary, depending on the size of the target quantity of interest (0D, 1D and 2D). Eventually, the deployment of the developed models in the context of an optimization problem will be tested.

In the following, the experimental test case used as a reference point for the present investigation is presented in Section II. The numerical setup is introduced and the ability in reproducing the results in agreement with the experimental data on the target configuration is demonstrated. Then, the design of experiments is introduced (Section III) and the methodology is deployed to generate the data (Section IV). The data are then used to train different algorithms depending on the fidelity level of interest: global quantities (0D), wall quantities (1D) or field quantities (2D). Sensitivity analyses are carried out to establish the best set of hyperparameters and loss function. Finally, an example of application of the developed models in the framework of an optimisation problem is shown in Section VI.

II. Reference test case: coaxial single injector

The database generated and analysed in the present work is based on a reference experimental test case [START_REF] Celano | Injector Characterization for a Gaseous Oxygen-Methane Single Element Combustion Chamber[END_REF]. It corresponds to a single coaxial injector chamber, operated at 20 bar, with gaseous methane and gaseous oxygen. The gaseous oxygen enters the square chamber of side 12 mm through the inner tube of the injector that has a radius of d o = 2 mm. The gaseous methane flows through the annulus around the oxygen, that has a height of d f = 0.5 mm.

RANS simulations are carried out with the commercial software Ansys Fluent with the objective of retrieving the stationary condition. As a first step, a numerical study has been carried out in order to find a numerical configuration that reproduces the experimental data at the reference point. This includes a mesh independence study, turbulence model sensitivity analysis and a comparison between several chemical schemes and turbulent combustion interaction models. The test case is largely documented and has been used as a benchmarck for several numerical studies, including RANS simulations [START_REF] Maestro | Large Eddy Simulation of Combustion and Heat Transfer in a Single Element GCH4/GOx Rocket Combustor[END_REF][START_REF] Chemnitz | Numerical Investigation of Reacting Flow in a Methane Rocket Combustor: Turbulence Modeling[END_REF]. In the following sections are given the details about the governing equations, the numerical setup and grid selected to carry out the present study.

A. Governing equations

The compressible, reactive, turbulent flow is described with the following conservation equations for mass, momentum and energy:

∂ ρ ∂t + ∂ ρ u j ∂ x j = 0 , (1)
∂ ρ u i ∂t + ∂ ρ u i u j ∂ x j = - ∂p ∂ x j + ∂τ i j ∂ x j , (2)
∂ ρ h ∂t + ∂ ρ u i h ∂ x i = ∂ ∂ x i λ + λ t c p ∂ h ∂ x i , (3)
where ρ is the density, p is the pressure, u i the velocity field component along the x i axis, λ and λ t are laminar and turbulent conductivities and c p is the specific heat at constant pressure. h is the specific total enthalpy (in J/kg) which is given by:

h = ũi ũi 2 + j ∫ T T r e f c p, j dT + h 0 j (T re f , j ) (4)
where h 0 j is the formation enthalphy of species j at the reference temperature T re f , j . The Favre formalism has been used and for every flow variable f , f is an ensemble average and f is a mass weighted average [START_REF] Pope | Turbulent Flow[END_REF][START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]. τ i j is the stress tensor which is given by the following expression:

τ i j = (µ + µ t ) ∂ u j ∂ x i + ∂ u i ∂ x j - 2 3 δ i j ∂ u k ∂ x k , ( 5
)
where µ is the laminar dynamic viscosity of the mixture and µ t is the turbulent viscosity. Thermophysical models are needed in order to close the system. The perfect gas equation is considered. NASA polynomials [START_REF] Mcbride | Computer Program for Calculation of Complex Chemical Equilibrium Compositions and Applications: 1. Analysis[END_REF] are employed for the specific heat at constant pressure c p with an ideal gas mixture law. Laminar conductivity and dynamic viscosity for the mixture are considered constant and equal respectively to λ = 0.0454W m -1 K -1 and µ = 17.2 µPa s. For what concerns turbulence closure modelling, the kmodel [START_REF] Launder | Lectures in Mathematical Models of Turbulence[END_REF] is used and therefore transport equations for the turbulent kinetic energy k and its rate of dissipation have to be solved in addition to the system of Eqs. ( 1)-( 3):

∂ ρk ∂t + ∂ ρ u i k ∂ x i = ∂ ∂ x i µ + µ t σ k ∂k ∂ x i + P k + ρ 2 k a 2 -1 , (6)
∂ ρ ∂t + ∂ ρ u i ∂ x i = ∂ ∂ x i µ + µ t σ ∂ ∂ x i + C 1 k P k -C 2 ρ 2 k , ( 7
)
where a is the sound speed. The production term of turbulent kinetic energy P k is expressed as follow:

P k = µ t ∂ u j ∂ x i + ∂ u i ∂ x j ∂ u i ∂ x j . (8)
The model constants used are C 1 =1.44, C 2 =1.92, σ k = 0.9 and σ = 0.9. The turbulent dynamic viscosity can be computed from k and :

µ t = ρC µ k 2 , (9)
where C µ = 0.09 is a constant. Eventually the turbulent conductivity k t is given by:

λ t = c p ρ µ t Pr t , (10)
where Pr t = 0.9 is the turbulent Prandtl number. The non-adiabatic Steady Diffusion Flamelet (SDF) approach is used to model the multi-species reactive flow [START_REF] Peters | Laminar Diffusion Flamelet Models in Non-Premixed Turbulent Combustion[END_REF], making use of a β-PDF to account for turbulence-chemistry interaction.

The Flamelet ). The GRI-3.0 mechanism, which was originally built for pressures around 1 atm, has been successfully applied to simulate higher pressure combustion in liquid rocket engines [START_REF] Chemnitz | Numerical Investigation of Reacting Flow in a Methane Rocket Combustor: Turbulence Modeling[END_REF][START_REF] Kim | Effects of Pressure and Inlet Temperature on Coaxial Gaseous Methane/Liquid Oxygen Turbulent Jet Flame under Transcritical Conditions[END_REF]. Therefore, it is considered suitable for the present investigation. A mixture fraction Z is introduced:

Z = Y i -Y i,o Y i, f -Y i,o (11)
where Y i is the mass fraction for the i th species, and the subscript f and o refer to the fuel and oxydizer inlets. Density ρ, temperature T and species mass fractions y i (generically called φ in the following ) are estimated from the local average mixture fraction Z, its variance Z 2 and the local average total enthalpy h (described by (Eq.3)):

φ = ∫ 1 0 φ( Z, h)P( Z, Z 2 )dZ , (12)
where P( Z, Z 2 ) is a standard β-PDF function. Therefore, two additional transport equations for Z and Z 2 are considered:

∂ ρ Z ∂t + ∂ ρ u i Z ∂ x i = ∂ ∂ x i ρ(D + D t ) ∂ Z ∂ x i , (13)
∂ ρ Z 2 ∂t + ∂ ρ u i Z 2 ∂ x i = ∂ ∂ x i ρ(D + D t ) ∂ Z 2 ∂ x i + C g µ t . ∂ Z ∂ x i 2 -C d ρ k Z 2 , ( 14
)
where C g = 2.86 and C d = 2. D and D t are laminar and turbulent diffusion coefficients respectively. Considering a unitary number of Lewis, Le=1, the laminar diffusion is expressed as a function of λ and µ:

D = λ ρc p . (15)
A turbulent Schmidt number Sc t = 0.5 is introduced in order to express D t as a function of µ t according to:

D t = µ t ρSc t . (16)

B. Numerical setup

To solve the above governing equations with the Ansys Fluent code, a pressure based implicit scheme, accounting for compressibility effects, has been selected to carry out axisymmetric simulations looking for a steady solution. The experimental chamber is rectangular and the nozzle has a rectangular cross section. In order to carry out axisymmetric simulations, the geometry has thus been modified into a cylindrical chamber and a conical nozzle. The radius of the chamber is d c = 6.77 mm , which corresponds to the same cross section area of the experimental square chamber (namely d c = A c /π with A c cross section area of the chamber). A schematic of the computational domain is reported in Fig. 1, with a close-up view of the injector and nozzle regions. Inlet boundary conditions are specified in terms of mass flow rate m and temperatures for both propellants, while no boundary condition is applied at the nozzle outlet where the flow is expected to be supersonic. Specifically, oxygen and methane are at 278 K and 269 K during injection respectively, the mass flow rate is m = 0.062 kg/s and the oxidizer-to-fuel ratio is (O/F) = 2.62.This translates to inlet velocities of 122 m/s and 132 m/s for oxygen and methane respectively. The walls of the injector assembly are assumed adiabatic with a no slip boundary condition. A no slip, isothermal boundary condition is applied on the chamber wall, considering a thermal resistance for a solid wall of copper with a thickness of 1 mm. The outer temperature profile is a linearization of the experimental wall temperature provided in [START_REF] Celano | Injector Characterization for a Gaseous Oxygen-Methane Single Element Combustion Chamber[END_REF] in the same position, as reported in Fig. 2.

x [mm]

T w [K]

C. Grid convergence analysis

The mesh is not uniform and the resolution is increased in the region where stronger gradients are expected. Indeed, a refinement is applied around the area of fuel and oxidizer mixing as can be observed in the detailed view of Fig. 1a. A grid sensitivity study was carried out in order to verify the requirements in term of resolution. Fig. 3 shows the wall heat flux profiles obtained with three different grids: coarse (181×10 3 cells), medium (720×10 3 ) and fine (2 770×10 3 cells).

The ratio of the cells height from one grid to another is 2. Wall heat flux profiles are similar, with differences between coarse and fine meshes smaller than 8%. In Fig. 3 are also reported the experimental curve [29] ( ) and the numerical results provided by [START_REF] Chemnitz | Numerical Investigation of Reacting Flow in a Methane Rocket Combustor: Turbulence Modeling[END_REF] (--). The present numerical results are in agreement with experimental data and present the same level of errors obtained with other simulations carried out with similar models in [START_REF] Chemnitz | Numerical Investigation of Reacting Flow in a Methane Rocket Combustor: Turbulence Modeling[END_REF]. From a quantitative point of view, the error between the solution and an extrapolated exact solution can be evaluated from the results obtained with the three grids, using the Richardson method [START_REF] Roche | Verification of Codes and Calculations[END_REF]. As an example, the error for the integrated wall heat flux Q w , the average total pressure p c and temperature T c at the chamber outlet (nozzle inlet) versus the number of grid cells are reported in log scale in Fig. 4. The error decreases by increasing the number of cells, showing a convergence of the results.

x [mm] In conclusion the convergence analysis is satisfactory and the obtained results are in agreement with experimental data. Considering that the purpose of the present work is to put in place a methodology by generating a large database, it has been decided to use the coarse mesh (181 000 cells) in order to reduce the computational time for the data generation.

q w [MW/m 2 ] [30]
All the results presented in the following are therefore carried out with a grid similar to the coarse mesh, with variations according to the geometry parametrisation.

III. Design of experiments and data generation A. Inputs and outputs

The reference numerical setup detailed in Section II will be used to generate a large database to infer models for injector design optimization. To do so, the first task is to define the inputs and outputs. The outputs are the quantities needed to define the best injector design. Therefore, these outputs are related to both performances and the sizing of other rocket engine components that depend on the injector design. Moreover, a visualisation of the temperature field inside the chamber can be useful to understand how the design choice impacts the flame and flow topologies. All the outputs are directly deduced from the results of each simulation. 0D outputs, referred as global quantities in the following, are computed with volume averages. 1D quantities are the wall quantities. 2D quantities are field quantities.

More information about these outputs are given in the next sections.

Regarding the inputs, seven geometrical parameters are necessary to define the geometry of the single coaxial injector chamber. These parameters, shown in Figure 1 are: the chamber radius d c , the contraction ratio between the chamber (c) and the throat (t) areas A c /A t , the expansion ratio between exit (e) and throat areas A e /A t , the oxygen inlet radius d o , the tip width d t , the fuel inlet height d f and the recess length l r . The only geometrical parameter which is kept constant for all the simulations is the chamber length l c , which has been increased with respect to the experimental value and is equal to 580 mm. This has been done in order to ensure a combustion length higher than the flame length for all the simulations. Two more inputs are used in order to define the operating point: the total inlet mass flow rate m and the oxidizer to fuel mass ratio O/F. The propellants injection temperatures, the wall temperature and the outlet atmospheric pressure (p a = 1 atm) are kept constants for all the simulations. To summarise, the design of experiments (DOE) will be build varying 9 parameters: 7 geometrical parameters and 2 parameters defining the operating point.

B. High-dimensional sampling

When sampling the 9D parameter space, it is required that design points are well spaced and give a good representation of possible parameters combinations. However, a purely random sampling along each axis does not guarantee a maximized spacing between points. Therefore, Latin Hypercube Sampling (LHS) [START_REF] Mckay | A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code[END_REF] is used to sample the parameter space. This method takes into account all previously sampled points when determining the next one. Specifically, the Enhanced Stochastic Evolutionary [START_REF] Jin | An Efficient Algorithm for Constructing Optimal Design of Computer Experiments[END_REF] algorithm which is an efficient implementation adhering to the LHS optimally criteria, is used with the Surrogate Modeling Toolbox (SMT) [START_REF] Bouhlel | A Python Surrogate Modeling Framework with Derivatives[END_REF]. 25), [START_REF] Chemnitz | Numerical Investigation of Reacting Flow in a Methane Rocket Combustor: Turbulence Modeling[END_REF], [START_REF] Pope | Turbulent Flow[END_REF] and [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF] in Appendix VII.A). Note that lower inlet sizes and chamber diameters are more likely to result in violation of these parameter limits. If these parameter ranges are sampled in a linear space, this will result in few design points with a combination of small inlets and chamber sizes. To increase the likelihood of these edge cases being sampled, d o , d t and d f are sampled in √ -space and d c is sampled in 4 √ -space. This implies that samples are taken linearly from functions f , g image spaces such that f (x) = √ x and g(x) = 4 √ x, with x representing the sampled variable. This can be observed in Fig. 5 that shows the histograms of the sampling points for the different parameters. The parameters sampling ranges and the sampling mode are summarised in Tab. 1. 6 Probability Density Function isosurfaces for all parameter combinations. x and y-axis limits as in Table 1.

0 d o d t d f d c O/F ṁ l r A c /A t A t /A e d o d t d f d c l r < l a t e x i t s h a 1 _ b a s e 6 4 = " G z G o g n i L T 3 + 7 U Z 2 H x H S p q 4 c o y y k = " > A A A C z n i c j V H L S s N A F D 2 N r 1 p f V Z d u g k V w V R I R d V k U x J 0 V 7 A P a I s l 0 W o N p E i a T Q g n B r T / g V j 9 L / A P 9 C + + M K a h F d E K S M + e e c 2 f u v W 7 k e 7 G 0 r N e C M T e / s L h U X C 6 t r K 6 t b 5 Q 3 t 5 p x m A j G G y z 0 Q 9 F 2 n Z j 7 X s A b 0 p M + b 0 e C O y P X 5 y 3 3 7 k z F W 2 M u Y i 8 M r u U k 4 r 2 R M w y 8 g c c c S V S n O x A O S y + z 9 D y 7 K V e s q q W X O Q v s H F S Q r 3 p Y f k E X f Y R g S D A C R w B J 2 I e D m J 4 O b F i I i O s h J U 4 Q 8 n S c I 0 O J v A m p O C k c Y u / o O 6 R d J 2 c D 2 q u c s X Y z O s W n V 5 D T x B 5 5 Q t I J w u o 0 U 8 c T n V m x v + V O d U 5 1 t w n 9 3 T z X i F i J W 2 L / 8 k 2 V / / W p W i Q G O N E 1 e F R T p B l V H c u z J L o r 6 u b m l 6 o k Z Y i I U 7 h P c U G Y a e e 0 z 6 b 2 x L p 2 1 V t H x 9 + 0 U r F q z 3 J t g n d 1 S x q w / X O c s 6 B 5 U L W P q t b V Y a V 2 m o + 6 i B 3 s Y p / m e Y w a L l B H Q 3 f 8 E U 9 4 N u r G 2 M i M + 0 + p U c g 9 2 / i 2 j I c P i S 6 T 6 w = = < / l a t e x i t > O F < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 u K Q j 7 D l F j j z G X y k 8 0 k / a c l l Z C 0 = " > A A A C y n i c j V H L S s N A F D 2 N r 1 p f V Z d u g k V w V V I R d V l 0 4 8 J F B f u A t k i S T m t o k g m T i V B C d / 6 A W / 0 w 8 Q / 0 L 7 w z T k E t o h O S n D n 3 n D t z 7 / W S M E i l 4 7 w W r I X F p e W V 4 m p p b X 1 j c 6 u 8 v d N K e S Z 8 1 v R 5 y E X H c 1 M W B j F r y k C G r J M I 5 k Z e y N r e + E L F 2 / d M p A G P b + Q k Y f 3 I H c X B M P B d S V S 7 N + A y j 6 a 3 5 Y p T d f S y 5 0 H N g A r M a v D y C 3 o Y g M N H h g g M M S T h E C 5 S e r q o w U F C X B 8 5 c Y J Q o O M M U 5 T I m 5 G K k c I l d k z f E e 2 6 h o 1 p r 3 K m 2 u 3 T K S G 9 g p w 2 D s j D S S c I q 9 N s H c 9 0 Z s X + l j v X O d X d J v T 3 T K 6 I W I k 7 Y v / y z Z T / 9 a l a J I Y 4 0 z U E V F O i G V W d b 7 J k u i v q 5 v a X q i R l S I h T e E B x Q d j X z l m f b e 1 J d e 2 q t 6 6 O v 2 m l Y t X e N 9 o M 7 + q W N O D a z 3 H O g 9 Z R t X Z S d a 6 P K / V z M + o i 9 r C P Q 5 r n K e q 4 R A N N X e U j n v B s X V n C m l j 5 p 9 Q q G M 8 u v i 3 r 4 Q N u Y J J O < / l a t e x i t > ṁ < l a t e x i t s h a 1 _ b a s e 6 4 = " L p M I 0 T T Y 4 h f L D 4 C b 3 0 Y 8 F 0 N V N q g = " > A A A C 2 H i c j V H L S s N A F D 2 N r 1 p f 1 S 7 d B I v g q q Q i 6 r L q x m U F + 8 C 2 l G Q 6 r a F p E i Y T o Z S A O 3 H r D 7 j V L x L / Q P / C O 2 M K a h G d k J k z 5 9 5 z Z u 5 c J / T c S F r W a 8 a Y m 1 9 Y X M o u 5 1 Z W 1 9 Y 3 8 p t b 9 S i I B e M 1 F n i B a D p 2 x D 3 X 5 z X p S o 8 3 Q 8 H t k e P x h j M 8 U / H G D R e R G / i X c h z y z s g e + G 7 f Z b Y k q p s v t P v C Z p O T 7 o Q l i V p k k n T z R a t k 6 W H O g n I K i k h H N c i / o I 0 e A j D E G I H D h y T s w U Z E X w t l W A i J 6 2 B C n C D k 6 j h H g h x p Y 8 r i l G E T O 6 R 5 Q L t W y v q 0 V 5 6 R V j M 6 x a N f k N L E L m k C y h O E 1 W m m j s f a W b G / e U + 0 p 7 r b m F Y n 9 R o R K 3 F N 7 F + 6 a e Z / d a o W i T 6 O d Q 0 u 1 R R q R l X H U p d Y v 4 q 6 u f m l K k k O I X E K 9 y g u C D O t n L 6 z q T W R r l 2 9 r a 3 j b z p T s W r P 0 t w Y 7 + q W 1 O D y z 3 b O g v p + q X x Y s i 4 O i p X T t N V Z b G M H e 9 T P I 1 R w j i p q 5 D 3 G I 5 7 w b F w Z t 8 a d c f + Z a m R S T Q H f h v H w A a k Q l 9 4 = < / l a t e x i t > A c A t < l a t e x i t s h a 1 _ b a s e 6 4 = " 1 f e O k I + N / I v 1 4 l u a e T l n 2 5 Y J u v Y = " > A A A C 2 H i c j V H L S s N A F D 2 N r 1 p f 1 S 7 d B I v g q q Q i 6 r L q x m U F + 8 C 2 l G Q 6 r a F p E i Y T o Z S A O 3 H r D 7 j V L x L / Q P / C O 2 M K a h G d k J k z 5 9 5 z Z u 5 c J / T c S F r W a 8 a Y m 1 9 Y X M o u 5 1 Z W 1 9 Y 3 8 p t b 9 S i I B e M 1 F n i B a D p 2 x D 3 X 5 z X p S o 8 3 Q 8 H t k e P x h j M 8 U / H G D R e R G / i X c h z y z s g e + G 7 f Z b Y k q p s v t P v C Z p O T 7 o Q n i V p k k n T z R a t k 6 W H O g n I K i k h H N c i / o I 0 e A j D E G I H D h y T s w U Z E X w t l W A i J 6 2 B C n C D k 6 j h H g h x p Y 8 r i l G E T O 6 R 5 Q L t W y v q 0 V 5 6 R V j M 6 x a N f k N L E L m k C y h O E 1 W m m j s f a W b G / e U + 0 p 7 r b m F Y n 9 R o R K 3 F N 7 F + 6 a e Z / d a o W i T 6 O d Q 0 u 1 R R q R l X H U p d Y v 4 q 6 u f m l K k k O I X E K 9 y g u C D O t n L 6 z q T W R r l 2 9 r a 3 j b z p T s W r P 0 t w Y 7 + q W 1 O D y z 3 b O g v p + q X x Y s i 4 O i p X T t N V Z b G M H e 9 T P I 1 R w j i p q 5 D 3 G I 5 7 w b F w Z t 8 a d c f + Z a m R S T Q H f h v H w A a 3 i l + A = < / l a t e x i t > A e A t uniform
d o d t d f d c l r d t d f d c O/F A c /A t A e /A t l r PDF Fig.
To further investigate the distribution of the design points, Fig. 6 shows a 2D projection of the probability density function (PDF) for all parameters combinations. The influence of the parameters limits is clearly visible. For example, the third row from the top with d c on the y-axis illustrates the dependence of Eq.( 32) on d o , d t and d f . The distribution of d c is however also influenced by Eq. ( 25), which results in a relationship between the contraction ratio and the chamber diameter seen in row A c /A t column d c . The contraction ratio is also influenced by d o , d t and d f via the chamber pressure approximation (Eq. 28). Other parameters like O/F, m, l r and A t /A e do not feature prominently in the limit equations and therefore do not lead to significant changes in the uniformity of the distribution. This can be observed by looking at their respective rows which feature little to no change in the distribution along the y-axis. In general, if the distribution evolves only in x or y direction, the two parameters have fewer dependence on each other via parameter limits. The more the gradient features a combined x-y direction, the more the two parameters interact in the DOE. Note that this only highlights the correlations between input variables when sampling the database, but is not related to any later correlation in the learnt surrogate model.

IV. Data generation, overview and preparation for training

Following the previous rules, a DOE of 3640 points is considered and for each point a RANS simulation is carried out until convergence using 16 cores on Xeon-6126 processors at 2.6GHz. On average, a run takes 52 minutes (13.87 CPU hours ) for a total of 50 000 CPU hours for the whole dataset generation. Over the total 3640 simulations, outliers outside of the 3σ range are removed which result in N = 3437 simulations that remain for training. Prior the training phase, a post-processing is necessary to extract outputs of interest from the generated data: these consist in global, wall and field quantities. These are described and analysed in the following sections.

A. Global quantities

Twelve global quantities of interest for the injector design will be computed: the average total pressure at the chamber outlet p c , the ratio between the total pressure loss and the inlet total pressure in the fuel injector ∆p f and in the oxidizer injector ∆p o , the maximum average temperature along the chamber axis T max , the flame length L f l estimated as the position of T max , the maximum wall heat flux q w,max , the average wall heat flux q w , the velocity magnitude at the nozzle outlet V E , the thrust F, the thrust coefficient c F , the characteristic velocity c * and the specific impulse I sp .

An overview of the ranges and distributions of the outputs is given by the histograms of Fig. 7. The ranges of these outputs were not known in advance and are reasonable with respect to the parameters ranges. For instance, the means for the chamber pressure p c and the peak wall heat flux q w,max are p c = 31.82bar and q w,max = 7.6MW/m 2 respectively, which are well in the same order of magnitude of the reference case (p c ≈ 20bar and q w,max ≈ 5.5MW/m 2 ). uniform 0

p c [MPa] Δp 0 [%] Δp f [%]

B. Wall quantity

The only wall quantity considered is the chamber top wall heat flux distribution q w . The wall profile is extracted from the simulation and sampled at 1024 equidistant points along the wall, making an interpolation from the computational grid. For the 3437 available samples, this amounts to a training dataset size of 30MB, which is 2 orders of magnitude larger than the size for global quantities. The heat flux exhibits a peak in the first 5-30 cm of the chamber and then The flame topology is strongly affected by the parameters choice. The overall range of temperature is similar and varies between the inflow temperature around 250 K and a maximum temperature between 2600 K and 3300 K.

The field data provided by CFD on the computational non-uniform mesh must be processed into a uniform Cartesian grid with a pixel-like format prior to the training. In fact this is a requirement in order to use convolutional neural networks (CNN). The fields are therefore interpolated on a uniform fixed resolution grid of size 128×128. The fixed number of interpolation points results in different spatial resolutions in the transverse direction of the injector depending on the chamber diameter. In Fig. 9 are reported three examples of temperature fields from three different DOE samples.

In the left part, the fields on the grid from the CFD simulation are shown. Note that the transverse y axis has been magnified to allow a better visualisation, being the chamber length always much greater than the diameter. In the right part of the figure are shown the corresponding transformed views of the same fields once interpolated on the 128×128 grid. These are the snapshots that will be used to train the CNN. In addition, a Boolean mask of the geometry is provided at the same resolution of 128×128. Specifically, the Boolean mask is a 128×128 table that assigns true values only in cells containing the fluid field, and false values in cells outside the fluid domain. The error will be computed only in the masked subset of the chamber and therefore only the cells containing the temperature information will be used for training. Additionally, all 9 scalar DOE parameters are supplied by nine constant 128×128 matrices: for each matrix, all entries are equal to the respective DOE parameter.

T [K]

Both the 9 DOE parameters and the temperature field are scaled to 0 mean and 1 standard deviation.

V. Results

Depending on the dimension of the output expected from the surrogate models (0D global quantities, 1D wall quantities or 2D field quantities) different algorithms will be tested for the training. Data are split in samples for training (2/3) and samples for validation (1/3). Furthermore, 77 samples are held back and never used during training. They will be used as test data in order to prevent overfitting when optimizing the hyper-parameters to the validation data.

A. Global quantities

Global quantities values are the training targets with the lowest complexity in terms of overall dataset size and target dimensions. This translates to low computational cost when training models and the necessity to use simpler models to avoid the possibility of over-fitting. Training data always include all DOE parameters, but only one output target at the time: the models are trained on the 12 post-processing quantities described in Section V. In the following, we will indicate with ŷ ∈ R n y the prediction of a model, which is generically given by:

ŷ = f (x, x t , y t ) (17)
where x ∈ R n x is the vector of inputs, x t ∈ R n x ×N the vector of training inputs and y t ∈ R n y ×N the vector of training outputs, where n x , n y and N are respectively the dimensions of the input, output and the number of samples. Therefore, for every global quantity it is n x = 12 and n y = 1. Three different algorithms are tested and compared for the surrogate model f . The first two are classical machine learning methods: Radial Basis Functions (RBF) and Kriging (KRG). For RBF, the following interpolating function f is used, with Gaussian basis functions:

f = N i exp - x -x t i 2 2 H 2 w R , (18)
where the coefficients vector w R ∈ R N is evaluated during the training, and the scaling parameter is set to H = 0.05.

For KRG, Gaussian covariance functions are used together with a linear model for the deterministic linear regression term, thus leading to the following expression for f :

f = β 0 + β 1 x + G(x) , (19)
where

β 0 ∈ R and β 1 ∈ R n x . G(x)
is the realisation of a stochastic process with zero mean value and a spatial covariance given by a Gaussian function:

cov[G(x i , x j )] = σ 2 Π n x k=1 exp -θ k (x k i -x k j ) 2 ) ( 20
)
where σ 2 is the variance of the process G(x) and the n x hyper parameters θ k are determined solving an optimization problem in order to obtain a maximum likelihood function. They are all initialised to 0.01. Both RBF and KRG are provided by the Surrogate Modeling Toolbox (SMT) open source package, and the reader can refer to [START_REF] Bouhlel | A Python Surrogate Modeling Framework with Derivatives[END_REF] for the details about the algorithms.

In addition, a fully connected neural network (FCNN) is also trained using the PyTorch library [START_REF] Paszke | PyTorch: An Imperative Style, High-Performance Deep Learning Library[END_REF]. A standard approach for regression is used with a ReLU activation function on all layers but the last, which has no activation function [START_REF] Chollet | Deep Learning with Python[END_REF]. The Adam optimizer is used with the Mean Squared Error Loss (MSE), defined by:

MSE(y, ŷ) = 1 N N i=1 (y i -ŷi ) 2 (21)
where y i is the target output and ŷi is the network prediction. The quality of a network is evaluated considering its accuracy as a figure of merit, which is related to the mean relative error over the samples of the validation dataset.

Specifically, the relative error for a target variable y i , for sample i, is defined with the following expression:

ε i = |y i -ŷi | y i , (22)
and the accuracy η is defined accordingly by:

η = 1 - 1 N N i=1 ε i × 100. ( 23
)
In order to build the suitable FCNN, two tuning studies are carried out. The first one aims at choosing the network architecture. A set of 23 different architectures with depths ranging from 1-4 hidden layers with layer widths between 16-1024 were tested with chamber pressure as target output. Given that preliminary tests had shown a tendency for overfitting even in small networks, all 23 architectures were trained with a weight decay of 1 × 10 -4 and a dropout of 0.01. Batch size was set at 64 with total training time of 500 epochs. 10-fold cross validation was used to allow for comparability between runs. Training these setups with 10-fold cross validation took roughly 3.5 hours on a 4-core machine. Given these results, the network architecture with the best combination of accuracy and lowest parameter count was chosen. This network has 2 hidden layers with 64 neurons each, so the overall architecture is 9 -64 -64 -1 with 4865 trainable weights and biases. The results obtained with the selected FCNN are compared against KRG and RBF models. It is worth mentioning that the default hyper-parameters values were chosen for RBF and KRG. To achieve a high confidence in the accuracy values obtained during training, 20-fold cross validation is used and the accuracy η on the validation dataset is used as a figure of merit and reported in Table 3. The best model for each parameter is highlighted in bold. The RBF setup never produces the most accurate model, and the corresponding accuracy is usually within 1 -2 % lower than the KRG accuracy. Generally, FCNN and KRG accuracies are in the same order of magnitude.

Parameter

As an example, Fig. 10a shows the predictions for training and validation data of both the KRG and FCNN models side-by-side for L f l and the corresponding relative errors ε. A slightly higher mean accuracy of 97.53% is obtained with the FCNN versus a 95.55% with KRG on the validation dataset (see Table 3). This difference mainly presents itself in the low end of combustion length, with comparably low relative errors for higher values. Training data error for KRG is not shown in the plot, as the KRG error in the training data is 0. Figure 10b shows the predictions for training and validation data of both the KRG and FCNN models side-by-side for q w,max . The FCNN again narrowly beats KRG with 97.50% mean accuracy versus 96.64%. The relative error of the FCNN model is consistently lower over the entire range of the target dimension. Finally, KRG and FCNN are compared in terms of inference times. Fig. 11 shows the inference time for a single predicted value, if it is part of a batch inference. On CPU, FCNN is around 10 times faster than KRG. A plateau is reached for batch sizes around 16 for KRG and 256 for FCNN. The fastest inference can be achieved for FCNN using a GTX1060 GPU. Moreover, the same batch inference time is obtained no matter the batch size until 65536.

time [ms]

batch size

B. 1D quantities: wall heat flux

The FCNN approach was again used to predict the heat flux all along the chamber top wall. The best four networks characteristics are summarised in Table 5. Most of the differences in peak accuracy for an architecture is due to the higher numbers of trainable weights leading to higher accuracy, albeit necessitating the use of higher dropout to prevent over-fitting. The mean accuracy is computed considering an average of Eq.( 23) over the wall length. The mean accuracy for training, validation and test data shows consistent values of around 96%. However, this mean accuracy for the entire length of the combustion chamber is not necessarily meaningful if the error occurs in critical locations like peak heating and on sharp gradients along the chamber. Fig. 13 illustrates the mean error distribution across all datasets and along the chamber. There is a strong correlation of the error distributions between all datasets. In agreement with the samples shown in Fig. 12, the larger mean relative errors occur early at the chamber inlet. The error quickly decreases thereafter and reaches lowest values in the last quarter of the chamber. Note that the deep learning models could be trained specifically to reduce this error if needed, for example by including the error on gradients into the loss function. Such a strategy is proposed in Section V.C, where it is proved effective to capture the sharp gradients of the temperature fields induced by the thin flame front.

Table 6 shows the mean relative error for the peak wall heat flux q w,max and the average wall heat flux q w as calculated from the wall heat flux curve. They translate to a mean accuracy η on validation data of 97% and 98%

respectively. These physical quantities are a good descriptors of the quality of the surrogate model, as they quantify the accuracy on typical quantities of interest. Comparing them to the accuracies obtained by the global modelling approach (reported in Table 3), both values are between 0.5% and 1% worse than the dedicated global model.

The 1D wall heat flux model under scrutiny was trained with a MSE loss function without any pondering of such quantities as q w,max and q w . In spite of this, the derived quantities from the model 1D predictions show to be in good agreement with those of the evaluated datasets. Recovering these quantities without specifically accounting for them in the loss function suggests the model is capturing the essential structures of the wall heat flux.

C. 2D quantities: temperature field

The entire temperature field in the injector and combustion chamber is learnt to provide in real-time a view of the complete field. A 2D field will be easier to interpret in terms of flame anchoring and shape, compared with 0D and 1D quantities. Therefore the output is a 128 × 128 = 16384 cells field. This 2D regression task was tackled with two variants of Convolutional Neural Networks (CNN). The first one is a simple encoder-decoder network using only convolutional layers. The second is a U-Net [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF] architecture including batch normalization, up-sampling layers and skip connections. The use of CNN on 2D spacial data has been studied extensively in the context of image processing and has recently been extended to be used on simulation tasks [START_REF] Chen | Numerical Investigation of Minimum Drag Profiles in Laminar Flow Using Deep Learning Surrogates[END_REF][START_REF] Ajuria | Towards an Hybrid Computational Strategy Based on Deep Learning for Incompressible Flows[END_REF][START_REF] Alguacil | Predicting the Propagation of Acoustic Waves Using Deep Convolutional Neural Networks[END_REF][START_REF] Thuerey | Deep Learning Methods for Reynolds-Averaged Navier-Stokes Simulations of Airfoil Flows[END_REF]. The setup of both architectures are displayed in Table 7. The channel and kernel sizing was partly based on networks in Thuerey et al. [START_REF] Thuerey | Deep Learning Methods for Reynolds-Averaged Navier-Stokes Simulations of Airfoil Flows[END_REF] and downscaled to fewer channels. The number of trainable weights is in the same order of magnitude as for the wall quantity learning task using FCNN networks. Learning rate and decay were also selected based on previously used values for this type of network [START_REF] Thuerey | Deep Learning Methods for Reynolds-Averaged Navier-Stokes Simulations of Airfoil Flows[END_REF]. Both architectures were trained using the same dataset of 3200 training, 170 validation and 77 test fields.

The uncompressed dataset size is 2.3GB with 10 128x128 inputs and 1 target image per simulation. First, a classical MSE loss, defined by Eq.( 21), is used. Note that this MSE is not computed over the whole 128 × 128 cells, but only on a masked subset containing the temperature information. The mask of the geometry provided as input is used for this purpose. The loss evolution over the number of epochs for both networks is shown in Figure 14. They both reach a plateau with regards to their validation loss after roughly 400 epochs of training. The U-Net training loss continues to decrease further, indicating an overfitting for larger epochs. This overfitting of the U-Net might come from its higher number of weights. No extensive hyper-parameters tuning was done for this task. samples. It can be seen that the amplitude of the error is less for U-Net. However, the location of the error remains similar. Both networks exhibit the higher relative error in very small areas of strong temperature gradients, associated with the thin flame front. The location of a strong gradient in temperature has a correlation with the combustion length, which is a quantity of interest for optimization. Therefore further action was taken to decrease the error in this region. A mix between MSE Loss and a Gradient Difference Loss (GDL) was investigated. The GDL loss is defined by:

GDL(y, ŷ) = 1 N N k=1 i, j || [y k (i, j) -y k (i -1, j)] -[ ŷk (i, j) -ŷk (i -1, j)] || 2 + ||(y k (i, j) -y k (i, j -1) -( ŷk (i, j) -ŷk (i, j -1)|| 2 . (24)
The GDL corresponds to an L2 loss on the gradient image, thus encouraging the network to capture correctly strong gradient regions. This strategy has been applied successfully in multiple fluid problems such as [START_REF] Alguacil | Predicting the Propagation of Acoustic Waves Using Deep Convolutional Neural Networks[END_REF]. Four different mixing ratios between the two losses were tested for both network architectures: only MSE, equal weighting (1/1) and increased weights for the GDL/L2 mix (10/1 and 100/1) . Figure 16 gives an overview of the relative errors. The left graph shows a decline in overall relative field error when using mixing ratios of 1/1 and 10/1 for both architectures, but a significant increase for 100/1. It suggests that an optimal mixing between GDL and L2 can be achieved. The optimal weighting 10/1 corresponds to a contribution of 79 % and 21 % of GDL and L2 in the total train loss, respectively.

Interestingly, the GDL allows for a better capture of the temperature gradients, i.e. a better accuracy on the flame shape. The field prediction can also be used to extract lower dimensional quantities known from the simulation data. Figure 17 shows the cross-section averaged temperature along the chamber for 4 different configurations with the same chamber length. It can be seen that the predicted curves closely match the true temperature distributions and most of errors occur at peak gradient values. The average temperature curves in Fig. 17 were evaluated to extract the peak temperature T max and its location, which is used as an indicator of the flame length L f l . These values were also obtained by the FCNN architectures trained on global quantities (see Sec. V.A). Figure 16 shows lower error for the U-Net architecture on these computed quantities. A mixed loss approach showed no improvement for T max . The mixed loss does however marginally reduce error for L f l . Indeed, a reduced gradient error might also yield to a reduced offset of gradient leading to better location accuracy down the chamber. Overall, the L f l error translated into a cell resolution is between 1 and 2 cells. A further decrease might therefore need higher resolution of the field.

p i n k c O 2 2 z f 2 T H i z v x B g 5 H 4 q w h l Y = " > A A A C 1 3 i c j V H L S s N A F D 2 N r 1 p f s S 7 d B I v g q i Q i 6 r L o x m W F v q Q t J U m n N T Q v k o m 0 h O J O 3 P o D b v W P x D / Q v / D O m I J a R C c k O X P u P W f m 3 m u F r h N z X X / N K Q u L S 8 s r + d X C 2 v r G 5 p a 6 X W z E Q R L Z r G 4 H b h C 1 L D N m r u O z O n e 4 y 1 p h x E z P c l n T G p 2 L e P O G R b E T + D U + C V n X M 4 e + M 3 B s k x P V U 4 u d g M J C n d a m v d Q z x 9 O e W t L L u l z a P D A y U E K 2 q o H 6 g g 7 6 C G A j g Q c G H 5 y w C x M x P W 0 Y 0 B E S 1 0 V K X E T I k X G G K Q q k T S i L U Y Z J 7 I i + Q 9 q 1 M 9 a n v f C M p d q m U 1 x 6 I 1 J q 2 C d N Q H k R Y X G a J u O J d
To summarise, the fact that using a GDL allows a better capture of the temperature gradients, is confirmed by the better accuracy obtained by the weighting 10/1 on both the whole temperature field as well as the location of the T max (that is L f l ). However, mixing GDL and L2 with 10/1 penalize the prediction of the maximal temperature in the combustion chamber. Moreover, this conclusion holds for both convolutional networks. It suggests that in general for multi-objective tasks, no proper loss can be defined, and a careful choice should be made depending on which objective is the priority. For instance here, a different mixing weight will be selected if T max or L f l is chosen as a priority. Finally, a field accuracy is estimated by computing the mean of Eq.( 23) over all the mask points for every sample and averaging over all the samples of the validation dataset. The accuracy is reported in Table 8 for Fully Convolutional and U-Net networks using a 10/1 GLD/L2 loss mix. The achieved accuracy with both networks exceeds 98%. Performances for T max and L f l are also reported. The accuracy on computed global quantities from the predicted field is comparable with the accuracy of the dedicated FCNN model. The U-Net architecture does have a slight edge over the Fully Convolutional network in all cases, the difference being most notable for L f l . While being more difficult and costly to train with respect to the FCNN on the global quantities, convolutional networks on the field data have the benefit to output the complete temperature field, much easier to interpret in terms of physics. Therefore, they constitute a first step towards reliable AI-based tools. error is also reported in order to highlight the sparsity of the results. For both 0D, 1D and 2D quantities, both the mean error and the sparsity of the error decrease when the number of samples N increases. However, the gain in increasing the DOE size is reduced for higher N. In particular, for N greater than 3000, the error for the global quantity reaches a plateau. This sensitivity analysis gives an idea of the requirements in term of DOE extension in order to obtain a specific accuracy. For instance, if an error between 1.3 and 2.3 % on the global quantities is acceptable, a dataset with N = 1000 samples is sufficient. On the other hand, if the surrogate model should provide high accuracy on the 1D and 2D quantities, the dataset used in the present analysis comprising N = 3437 samples, seems to be necessary. This analysis clearly points out the necessity in providing innovative strategies in order to reduce the dataset dimension if the methodology is to be applied on higher fidelity computations (namely Large Eddy Simulations). This will be the object of future works.

P

VI. Surrogate models application and exploitation

The surrogate models derived in section V can be used in the design and optimization of a coaxial injector and of the rocket configuration in which it is used. Specifically, these models can be useful in two ways. First, they can be used to carry out an optimization of the injector with respect to some specific figures of merit, and the optimization of a whole rocket system using the injector. In particular global and wall quantities are the more suitable outputs for this task.

Second, they can contribute to gain insight into some specific flow and combustion features that eventually lead to the specific performances that bring to the optimal design. Indeed, having access to the field with a very low inference time allows the immediate observation of the flow topology induced by a selected design and thus to deduce specific conclusion on the observed behaviour.

To showcase the items mentioned above, a constrained multi-objective injector design problem has been addressed.

It is of common interest for design bureaus to search for injector geometries which provide a best compromise between engine performance, typically represented by its specific-impulse I sp , and engine compactness. Both variables have an opposing dynamics as larger engines will typically allow for better combustion as flow residence time in the combustion chamber is increased. The dimensioning of the combustion-chamber length is associated to the characteristic length which directly depends on the flame length L f l [START_REF] Khan | Optimum characteristic length of gas generator for liquid propellant rocket engine[END_REF]. The rationale is smaller flame-lengths allow greater system compactness.

Fixed Values

Variables Ranges

F t [N] δ F [N] x = d o , d t , d f , d c , A c A t , A e A t T m [kg/s] O/F l r [mm]
135.9 1.0 [2.0, 0.5, 0.5, 6.77, 2.5, 2.0] 0.001 -0.13 2.0 -3.99 0.0 -30.0 The design problem under consideration is to find then x ∈ X p , where X p corresponds to the set of non-dominated solutions, also called Pareto solutions, of the multi-objective optimization (MOO) problem with objective functions:

f 1 (x) = -I sp (x, x
) and f 2 (x) = L f l (x, x). L f l corresponds to the axial location of the maximum average temperature in the chamber, hereby used as an indicator of flame-length. The problem is also subjected to the constraint g(x) ≤ 0, such

that g(x) = |F t -F(x, x)| -δ F .
Several algorithms which treat the MOO problem described above are available in the literature. In the current work, the NSGA-II (Non Dominated Sorting Genetic Algorithm) is employed. NSGA-II, as with other algorithms of its type, is non-gradient based. Note that even though NSGA-II spares of the readily available gradients offered by a neural-networks, the fast inference times of the surrogates models used (F, I sp , L f l ) allow to keep the optimization problem tractable. The algorithm is chosen due to its straight-forward implementation and accessibility. In the results to follow, the concerned values taken are detailed in Tab. 9.

The NSGA-II algorithm is run by means of the PyMOO python package [START_REF] Blank | Pymoo: Multi-Objective Optimization in Python[END_REF], with a population size of 1000 agents, random sampling and crossover probability of 0.9. The algorithm reached a solution set after 45 generations with a tolerance on the objective function of 10 -3 . Overall, 45000 evaluations where made of the objective functions f i (x).

The algorithm is observed to explore the complete range of the individual coordinates of the search-space. This is a convenient property of the algorithm which makes it more robust in the face of local minima, even though it is not a sufficient condition in the path towards the global minima. 9 which may result as a consequence of the strong relationship between the expected thrust and mass flow rate.

In the l r dimension the points are not distributed uniformly but rather three clusters are observed, as indicated in Fig. 19.

Cluster A corresponding to solutions where l r ≤ 2.0 mm, Cluster B where 20 mm ≤ l r ≤ 29.9 mm and Cluster C with l r > 29.9 mm. The non-dominated solutions labeled by cluster are plotted in the objective space in Fig. 20, also known as the Pareto front.

Points in Fig. 20 The errors for the I sp are 0.2%, 0.43 % and 0.52% respectively. These values are in the expected ranges discussed in Section V.A, since an average error of 3.47 % was reported for the L f l and 0.32% for the I sp . The impact of this error on the Pareto front estimation can be seen in Fig. 19 where the computed values of L f l and I sp for the three points from the RANS are also reported: points A, B and C from the RANS are in proximity of the Pareto front as expected from the estimated errors of the surrogate models.

The previous analysis showcased a particular application problem for multi-objective optimization with the aforementioned surrogate models, for a problem involving three design parameters. However, more intricate single objective optimization problems, involving a compromise objective function may be solved by harnessing the power of neural networks and the readily available gradients. Optimization problems, involving a larger number of design parameters could be addressed in future works. Combined with robust gradient based optimization algorithms such as Stochastic Gradient Descent and the fast inference time these models offer, multidisciplinary optimization loops may be solved in realistic timescales.

VII. Conclusions

In the present paper, a first attempt to jointly use data from computational fluid dynamics simulations and deep learning for the design of rocket engines injector and combustion chamber has been carried out. Reynolds-Averaged Navier-Stokes simulations of a single coaxial injector chamber, operated with gaseous methane and oxygen propellants, It is worth mentioning that from a quantitative point of view the presented results are of course affected by the validity and quality of the data. Reynolds-Averaged Navier-Stokes simulations do not allow to capture specific low scale flow rate can be estimated with the well known expression derived from the theory of 1D isentropic flow in nozzles [START_REF] Sutton | Rocket Propulsion Elements[END_REF]:

m = Γ p c A t √ RT c (26)
where Γ = γ 2 γ+1 γ+1 γ-1 , R is the gas constant and T c is the total temperature in the chamber. This equation can be recast in order to express the dependence of p c over three design parameters m, A c /A t and d c :

With this approximate pressure p c,ap , we are able to express the following constraints on d o , d f and d c in order to verify conditions 1) and 3):

d o > m o 2 0.7p c,ap π c v RT o 2 c p W O 2 , ( 30
)
d f > d 2 o + 2d o d t + d 2 t + m CH 4 0.7p c,ap π c v RT CH 4 c p W CH 4 -d o -d t , (31)
d c > a c a in d 2 o + d f 2d o + 2d t + d f . (32)
where R = 8314J/kg/K is the universal gas constant and W i is the molar mass of species i.

2 )

 2 Urbano § ISAE-SUPAERO, University of Toulouse, France Facing the need to increase the accuracy of rocket engines design tools, the present work introduces an innovative methodology for the design and optimisation of rocket engine combustion chambers using numerical simulations and deep learning. An experimental test case of a single coaxial injector is taken as a reference point and a design of experiments is generated by varying nine parameters (geometrical and operative conditions). Reynolds-Averaged Navier-Stokes simulations are carried out to generate the dataset. The data are used to train surrogate models of different fidelity, from low dimensional outputs (0D and 1D) towards the 2D temperature field. Attention is given on the selection of the proper machine learning technique. For low dimensional outputs, results show that deep neural networks outperform other standard machine learning tools, namely Radial Basis Function and Kriging. Regarding high-dimensional outputs, convolutional neural networks with gradient-based loss functions are found effective to capture the large and smooth temperature variations, as well as the thin and sharp temperature gradients at the flame front. Eventually the models are used in the framework of an optimisation problem. Results highlight the benefits of new design and optimisation tools based on deep learning, capable of real-time predictions of complex component along the x i axis [m/s] λ = thermal conductivity [W/m/K] a = speed of sound [m/s] c p = specific heat at constant pressure [J/kg/K] h = specific total enthalpy [J/kg] h 0 j = formation enthalpy of species j [J/kg] γ = specific heat ratio Γ = function of γ τ i j = fluid shear stress tensor [Pa] µ = dynamic viscosity [Pa.s] k = turbulent kinetic energy [J/kg] = kinetic energy dissipation rate [J/kg/s] P k = production of turbulent kinetic energy [= β-PDF function D = mass diffusion coefficient [m 2 /s] d = radius or height [m] l = length [m] A = cross section area [m 2 ] m = mass flow rate [kg/s] O/F = oxydizer to fuel mass ratio q w = wall heat flux [W/m 2 ] Q w = integrated wall heat flux [W] ∆p = relative total pressure drop [%] L f l = flame length [m] |V E | = velocity magnitude at the nozzle exit [m/s] I sp = specific impulse [s] F = thrust [N] c F = thurst coefficient c * = characteristic velocity [m/s] x = vector of inputs, coordinate y = vector of outputs N = number of samples in a dataset n x , n y = number of inputs and outputs ε = mean relative error of a model η = surrogate model accuracy w R , H = radial basis function parameters and coefficients β 0 , β 1 , θ K = Kriging model parameters δ F = uncertainty on the target thrust F t f , g, G = generic functions W = molar mass [kg/mol]

  Fuel and oxidizer inlet detail with recessed post tip.

  Refined nozzle and outlet detail.

Fig. 1

 1 Fig. 1 Parametrized geometry and meshes close to the injector tip (a) and nozzle (b). Note that the recess length is zero in the reference experiment ( l r = 0).

Fig. 2

 2 Fig. 2 Experimental wall temperature and linearization imposed as a boundary condition in all the simulations in the outer side of the 1mm thick chamber wall.

Fig. 3 Fig. 4

 34 Fig. 3 Wall heat flux profile. Comparison between experimental data [29], numerical results from [30], and the present setup for different meshes.

Fig. 5

 5 Fig. 5 Histograms (n Bin = 40) for all the parameters of the Design of Experiments (DOE) along their respective ranges: 1D projections of the 9D parameter space sampling.

Fig. 7

 7 Fig. 7 Histogram of the global quantities for N=3437 simulations. Outliers larger than 3σ have been removed.

  remains almost constant. The distribution of the heat fluxes over the whole DOE is shown in Fig.8. The minimum and maximum lines are combinations of different simulations results: these are the maximum and minimum heat flux at specific axial positions. The wall heat fluxes in rocket engines combustion chambers constitute a critical, but yet difficult to estimate, quantities during preliminary designs. Knowledge of these and their gradients is critical to the design of the chamber cooling system to ensure safe engine operation.

Fig. 8

 8 Fig. 8 Wall heat fluxes distributions over the Design of Experiments.

Fig. 9

 9 Fig. 9 Temperature field examples. Up: geometric view (y-axis magnified by a factor of 10). Bottom: transformed view used during training.

  Prediction (top) and relative error ε (bottom) versus true value for the flame length L f l . Prediction (top) and relative error ε (bottom) versus true value for the maximum wall heat flux q w, ma x .

Fig. 10

 10 Fig. 10 Top figures: model predictions for L f l and q w,max for training and validation data using KRG and FCNN. Bottom figures: relative error and moving mean.

Fig. 11

 11 Fig. 11 Inference times for a batch divided by batch size with trained KRG and FCNN models. Mean of 10 runs.

Fig. 12

 12 Fig. 12 Wall heat flux predictions (top) and corresponding relative errors (bottom) for five samples of the test dataset.

Figure 12 (

 12 Figure 12 (top) shows predictions for five randomly chosen samples from the test dataset. These samples have never been part of the training process of the network and are a good representation of the diversity of the dataset. The predictions were obtained using the network with the highest validation accuracy in Table5. The mean relative accuracy of 96.85% on the entire test dataset is comparable with the accuracy on the train and validation datasets, leading to the assumption that the network is able to generalize and has not overfit to training data or even to validation data during hyper-parameters tuning. The samples predictions show good accuracy on the overall shape of the wall heat flux curve along the chamber wall, as well as in capturing the sharp gradients and peaks in the first part of the chamber. The

Fig. 13 Mean

 13 Fig. 13 Mean Relative Error on the wall heat flux prediction along the chamber axis, for the three dataset (training, validation and test).

Fig. 14

 14 Fig. 14 Loss progression during training for base U-Net and Fully Convolutional (FC) networks , for training and validation datasets.

  An overview of the relative error distribution a test dataset prediction using either the fully convolutional net or the U-net is shown in Figure15(c) and (d), and the corresponding target output is shown in Figure15 (a). These figures illustrate the comparison only on one test dataset sample, taken as an example, but the behaviour is similar for all the Predicted temperature field using U-Net, 10/1 GDL . Relative error ε T using the Fully Convolutional Net. Relative error ε T using the U-Net. Relative error ε T using the U-Net with a mixed 10/1 GDL.

Fig. 15

 15 Fig. 15 Target and predicted temperature fields via best performing network, and relative errors fields for different networks. Sample of test dataset

For

  all cases, training, validation and test errors are close indicating only minimal overfitting. On average, the U-Net exhibits a small advantage for all runs compared to the Fully Convolutional architecture. An example of field prediction and corresponding relative error achieved with the best performing U-Net architecture using 10/1 mixed GDL/L2 Loss during training are shown in Fig.15(b) and (e) respectively. The strongest temperature gradient is well resolved. The relative error is concentrated in a very limited area in the fuel injector and as noise in the oxidizer injector lines. t e x i t s h a 1 _ b a s e 6 4 = "

TFig. 16

 16 Fig. 16 Mean relative error for trained CNN and U-Net networks with varying levels of GDL/MSE mixing in the loss. Top: T field. Middle: T max . Bottom: L f l .

Fig. 17

 17 Fig. 17 Mean temperature prediction along the chamber length. Predictions via U-Net (10/1 GDL).

Fig. 18

 18 Fig. 18 Error on the test dataset for different DOE sizes. Whiskers: complete range. Box: 1 st till the 3 r d quartile. Median line: mean error.

  Fig. 19 NSGA-II non-dominated solution set (Pareto solutions) frequency in the search-space.

Fig. 21

 21 Fig. 20 Projection of the non-dominated solutions in the objective space (Pareto front), from the multi-objective problem. Crosses: RANS results.

  have been carried out in order to generate a large database which comprises more than 3600 different points. It has been shown that particular attention must be devoted to the design of experiments generation in order to ensure the physical correctness of the computed configurations. The core of the work has been devoted to the development of three classes of surrogate models: global quantities, wall quantities and field quantities surrogate models. All the models have the same nine inputs parameters that define the geometry and the operative conditions and differ in the dimension of the output. A comparison between fully connected neural networks (FCNN) and other machine learning approaches, namely Kriging and radial basis functions, has been done for the gloabl quantities models, showing a comparable accuracy between the approaches. After hyperparameter tuning, the best neural network accuracy for the global quantities varies between 97.5 and 99.92%. When increasing the dimension of the problem (1D and 2D) only deep learning methods are of interest. For the wall heat flux, the developed surrogate model based on a fully connected neural network is able to reproduce the data with an accuracy up to 97%. Two types of convolutional neural network have been tested for the temperature field models: a standard CNN and a U-Net. Particularly, an increase in the accuracy is achieved by considering a mixed loss which combines a Mean Squared Error loss and Gradient Difference loss, that allows to reach a mean accuracy of 98.49% with the U-Net. Interestingly, the field surrogate model allows to retrieve important global quantities, like the flame length, with the same accuracy of the best FCNN selected for the corresponding global model. Eventually, the obtained models have been used in the framework of an optimisation problem. In this example, the global quantities models allow to choose the best configuration with respect to the objective function and the specific constraints. Immediately, the field model allows to visualise the corresponding temperature field thus giving clear indications and explanations of the observed behaviour in terms of performance.

2 ( 27 )

 227 /A t ) (A c /A t ) re f d c,re f d cwhere the index re f refers to the reference setup and K = ( Γ/RT c )/( Γ/RT c ) re f . From Eq.(27) we expect, in a first approximation, to find a linear relation between p c and m, and between p c and A c /A t and a quadratic relation between p c and d c .

Fig. 22 Fig. 22 .

 2222 Fig. 22 Response of chamber pressure p c to changes in parameters, sweeping through one parameter range at a time.

Table 1 Parameter range, reference values and sampling mode for all 9 parameters.

 1

	Parameter	Symbol	Range	Reference	Sampling Mode
	Oxidizer inlet radius	d o	0.5 -4 mm	2 mm	linear in	√	-space
	Injector tip width	d t	0.2 -1 mm	0.5 mm	linear in	√	-space
	Fuel inlet height	d f	0.2 -1 mm	0.5 mm	linear in	√	-space
	Chamber radius	d c	4 -10 mm	6.77 mm	linear in 4 √	-space
	Mixture ratio	O/F	2 -4	2.62	linear
	Total mass flow rate	m	0.005 -0.13 kg/s 0.062 kg/s	linear
	Recess length	l r	0 -30 mm	0 mm	linear
	Contraction ratio	A c /A t	2 -5	2.5	linear
	Expansion ratio	A t /A e	2 -3	2	linear

Some combinations of input parameters could lead to non-physical configurations (for instance a supersonic flow in the injector). To avoid non physical behaviours, some limits are imposed in terms of minimum of A c /A t , d o , d c and d f (see Eqs. (

Table 2 Hyper-parameters, with range and discretization, used during tuning for the FCNN (9-64-64-1) on the global quantities.

 2

		Learning Rate LR Decay Batch Size Weight Decay Dropout
	Range	0.001 -0.05	0.5 -0.99	8 -128	0 -0.002	0 -0.05
	Discretization	Uniform	0.01	8	0.0001	0.001
	Best setup	0.01214	0.97	32	0.0002	0.0
	The second analysis aims at finding the best hyper-parameters setup. The five hyper-parameters under consideration
	are the learning rate (LR), LR decay, batch size, weight decay and dropout. A Bayes search optimization approach
	is applied using the wandb sweep API. Three-fold cross validation is used leading to a 2/3-1/3 split between training
	and validation data per training run. All parameters ranges for tuning are listed in Table 2. Training would run for a

maximum of 200 epochs or if the learning rate hits 10 -10 . In total, 193 different setups were tested and the optimizer set to minimize the relative error on the validation dataset. Eventually, the best setup is summarized in Table

2

.

Table 3 Accuracy η on validation dataset (Eq.(23)) for global quantities models, obtained using 20-fold cross validation. Highest accuracies are highlighted in bold.

 3 For the three algorithms, training times are negligible compared to computational cost to generate the data.

	Specifically, training on a 4-core CPU took 1s for RBF, 52 s for KRG and 85 s for FCNN (considering 200 epochs and a
	32 batch size).

Table 4 Tested architectures and hyper-parameters ranges for training on wall heat flux.

 4

	Parameter	Layers	Learning Rate LR Decay Batch Size Decay Dropout
		9 -256 -256 -1024
		9 -64 -64 -128 -1024
	Range	9 -64 -128 -256 -1024	0 -0.1	0.5 -0.999	16 -256	0 -0.02	0 -0.1
		9 -256 -64 -256 -1024
		9 -128 -128 -128 -1024
	Discretization	Categorical	Uniform	0.001	16	0.0001	0.001

The main difference being that the last layer of the network now has 1024 nodes instead of just a single one. Due to the high dimensionality of the output data and the overall data size, KRG and RBF models were no longer pursued. Given the increased network size and gained experience from global model training, cross validation was no longer deployed. The dataset is split 90/10 between training and validation samples with again 77 samples held back as an additional test dataset. A tuning is carried out following the approach used in section V.A. In total, five different architectures are tested. The ranges of parameters are increased compared to the one used for the global quantities, as reported in Table

4

. Additionally, the maximum number of training epochs is increased to 1000.

Table 5 Best training setups during hyper-parameters tuning for wall heat flux. Mean accuracy evaluated for the training, validation and test datasests.

 5

Table 6 Mean Relative Error and Standard Deviation for peak

 6

(q w,max ) and average heat flux (q w ) across different datasets (training, validation and test).

Table 7 Convolutional Neural Network architectures used to predict temperature fields.

 7

	Type	Hidden Layers Trainable Weights Epochs Batch Size	LR	LR-Decay
	Fully Convolutional	9	400k	500	64	0.001	0.99
	U-Net	13	490k

Table 8 Test data accuracy for CNN architectures using 10/1 GDL mixed loss, compared against FCNN networks predicting global quantities.

 8

	P Quantity P P P	P	P Architecture P P P P P P	Fully Convolutional	U-Net	FCNN
	Temperature Field: T	98.35%	98.39%	-
	Peak y-Avg Temperature: T max	99.64%	99.73% 99.83%
	Location of T max : L f l	97.59%	97.85% 97.53%

Table 9 Multi-objective optimization problem variables values and ranges.

 9 The problem addressed is to find the subset of mass-flow rate m, mixture ratio (O/F) and recess length l r within the domain of validity of the global quantities models detailed in Section V.A which maximizes specific-impulse and minimize flame-length. Additionally, the design is constrained to yield a target thrust F t with a given uncertainty δ F .

Among the 9 input parameters of the DOE (see Sec III.A) only 3 are considered variables: m, O/F and l r . Remaining 6 degrees of freedom, hereby called x have been given arbitrarily fixed values (Tab. 9) to further simplify the present problem.

Funding Sources

This work was supported by the Chair for Space Advanced Concepts Laboratory (SaCLaB) resulting from the partnership between Airbus Defence and Space, Ariane Group and ISAE-SUPAERO. This work was supported by Région Occitanie, with the MODAR project (grant number: 20007323/ALDOCT-001030) and CERFACS for the financing of Jose Felix Zapata Usandivaras' PhD.

and non stationary phenomena that play an important role in the flow development inside the combustion chamber.

Moreover, all the study has been carried out considering a constant wall temperature, which is a strong approximation.

Nevertheless, the results show the potentiality of the methodology in permitting fast and comprehensive analysis of complex phenomena. This work represents a proof of concept of a methodology that needs to be further developed increasing the quality of the data. Specifically, it would be interesting to use Large Eddy Simulation to generate the data. Furthermore, it would be interesting to test the above developed surrogate models in the framework of a multi-disciplinary analysis and optimisation of a whole rocket engine.

Appendix A. DOE generation and constraints

Having chosen the input parameters, a database need to be generated by properly varying then inputs in the 9D parameter space, that is by generating a proper Design of Experiments (DOE). One possibility would be to randomly sample points linearly within the parameter space. This approach has three main drawbacks. First, some combinations of input parameters could lead to non-physical configurations and therefore have to be withdrawn from the DOE. Second, a random sampling would not guarantee a maximised space between the points. Third, a linear sampling could not be suitable depending on the mathematical dependence expected between inputs and outputs. Therefore, two improvements were made to increase the viability of the design points and their use in training a neural network as explained in the following. As expected, preliminary tests showed that a purely random sampling of the parameter space can lead to nonphysical setups. Mostly, these issues would manifest themselves as supersonic flows in the oxidizer or fuel inlets or chamber pressures orders of magnitude higher than anticipated. To avoid those points, we need to impose the following three conditions that are usual for the design of liquid rocket engine combustion chambers and injectors [START_REF] Sutton | Rocket Propulsion Elements[END_REF]. First, 1), the contraction area should be high enough in order to ensure a choked nozzle. Second, 2), the injector Mach number should not exceed 0.7: M in < 0.7. Third, 3), the chamber Mach number should not exceed 0.3: M c < 0.3 . Conditions 1) and 3) bring to the following restriction on the contraction area A c /A t :

with M c = 0.3, and γ is the specific heat ratio of the gas flowing in the nozzle. To find suitable constraints for d o , d f and d c in order to respect conditions 2) and 3) the pressure chamber p c influence has to be investigated. In fact, as the mass flow rate and temperature are inlet boundary conditions, the flow velocity in the inlet mainly depends on chamber pressure p c . The pressure chamber p c is an output of the simulation. In fact the nozzle is choked, and therefore the mass