Development of a robust solver to model the flow inside the engines for high-speed propulsion

The demand for discovering new commercial routes as well as the possibility to shortening civilian long-haul flights boosted the interest of civil hypersonic vehicle designs. Among all the multiple projects started by the various nations, the European community funded project STRATOFLY aims at refining the baseline LAPCAT II-MR2.4 design for further improvements. The new aircraft would enable a flight shorter that 3 hours from Brussels to Sydney, carrying 300-passengers above the already crowed atmosphere. The wide Mach range operability, up to Mach 8, demands the use of multiple engines, leading to a highly integrated propulsion system. The current study is focused on the development of new CFD platform to estimate the performance of the combined propulsion system during the supersonic to hypersonic transition. In order to control the complex flow physics, highfidelity CFD simulations remain the fundamental tools for the preliminary investigations. On the current framework, an advanced robust compressible solver has been developed in order to handle the different flow regimes. The new tool solves Unsteady Reynolds-Averaged Navier-Stokes (URANS) equations by employing cell-centered Finite Volume Method constructed on openFoam toolbox. Two innovative high-order discretization schemes, with different abilities, based on approximated Riemann solvers were developed for capturing the flow physics within high-speed propulsion systems. Advanced time discretization has been taken into account to increase the temporal accuracy. At the end, the whole implementation has been validated in multiple test cases, ranging from incompressible to hypersonic regimes, confirming its excellent stability, robustness and accuracy.


Introduction
Use During the last century, several international programs has started with focus on hypersonic vehicle designs. For instance, the Russian CIAM program performed the first test on a Dual Mode Ramjet (DMR), the Australian HyShot2 program [1] as well as the recent American Hyper-X program [2] demonstrated the possibility of a stable air-breathing hypersonic flight. The first European interest in civil supersonic and hypersonic flight was first represented by the European Community (EC) funded ATLLAS I project in 2008 [3]. Subsequentially, ATLLAS II [3] carried out the analysis concerning the integration of different propulsion system into a single supersonic vehicle. Progressively, LAPCAT I [4] project aimed of developing various vehicle concepts in order to potential reduce the flight time of antipodal destinations. The hypersonic concept was represented by the MR1 concept which, for the first time, was equipped by combined propulsion systems based on turbo cycle.
The last follow up of this project is represented by the LAPCAT II [4] which was the natural continuation of this concept. The LAPCAT II-MR1 was completely redesigned in order to increase the overall efficiency, including a more efficient combined propulsion system, leading to the LAPCAT MR2.4. At this point, the EC founded STRATOFLY project wants to refine the previous LAPCAT MR2. 4 concept with important investigations that could not be addressed within the timeframe of the previous projects. The idea is then to take as reference vehicle the LAPCAT MR2.4, making some changes in order to shorten the civil flight time of an order of magnitude with respect to the current state of art of civil aviation. Such investigations include several research areas, e.g. supersonic combustion, pollution, noise emission, lightweight structures and a deep investigation of the newer combined propulsion system. Therefore, the ability to fully control the complex flow physics of these advanced systems is fundamental to develop hypersonic vehicles. However, at this moment, no facilities are present to test full scale models and no appropriate scaling rules have been found. Moreover, realistic high-enthalpy conditions can only be obtained for a very short period [5]. High-fidelity CFD simulations remain then the fundamental tool to address the complex flow patterns. The implementation has to be as flexible as possible in order to be adaptable to every condition but accurate enough to avoid spurious results. The request for accurate CFD simulations is then mandatory to consolidate the final design. The aim of the project is to develop an open-source platform to simulate supersonic and hypersonic airbreathing propulsion systems. The choice of the most suitable schemes takes into account the accuracy, the computational cost and the versatility in order to build a robust solver for simulating various flow conditions. The originality of this work lies mainly into the combination of available numerical schemes to perform accurate highfidelity simulations of advanced propulsion systems in a modern computational platform. The implementation will be then verified on multiple test cases to establish its reliability, ranging from subsonic to supersonic. Subsequently, the tool will be used to evaluate the nozzle performances of the LAPCAT II-MR2.4 concept.

Numerical implementation
OpenFoam 4.0 was considered as reference platform for the implementation making use of the Finite Volume Method (FVM) thanks for its flexibility and high-accuracy for the problems involved. The cell-centered approach was then employed to solve the nonreactive Navier-Stoker equations by using conservative variables [6] in terms of densitybased architecture. Moreover, the assumption of calorically perfect gas was taken into account as well as the Fourier's law to discretize the diffusive energy flux [6]. The current version is not then able to solve real gas flow in which real thermodynamic process and chemical reactions take place. The bottom-line idea for the implementation is to merge the different peculiarities of available solvers as rhoCentralFoam, densityBasedTurbo, hyFoam [7] etc. in order to build a robust and reliable CFD tool to analyse high-speed air-breathing engines.

Time discretization
For the time discretization method of lines [6], which allow decoupling spatial and temporal discretizations was taken into account thanks to the larger flexibility. It basically reduces the governing equations to a system of coupled ordinary differential equation (ODE). At this early stage, the use of explicit time integration scheme was preferred since numerically cheap and requires only modest amount of memory. Moreover, to increase both stability and accuracy, a multistage forth-order Runge-Kutta time-stepping was employed. The idea is to advance the solution in few stages where the residuals are evaluated only on these different stages. The stability and the robustness are enforced by weighting the various coefficients and modifying the stage steps [6]. Moreover, such configuration is particularly suitable for upwind spatial discretization and it can be easily switched to a first-order scheme in case of strong shock. Finally, by proper tuning of the stage coefficients, this approach is appropriate for both steady and unsteady simulations.

Space discretization
The variation of the conservative variables on the numerical control volume depends only on the evaluation of inviscid and viscous fluxes through the surrounding surfaces. Concerning the viscous fluxes, due to the elliptic nature, central averaging was employed. More in details, the evaluation of the gradients was computed making use of the Green's theorem [6] performed by face averaging. On other hand, for the inviscid flux discretization, the use of upwind schemes as approximate Riemann solvers seemed to be efficient in terms of accuracy and computational cost. Indeed, upwind schemes are more advanced since they are constructed by considering the physical properties of the Euler equations and then a lower sensitivity to grid distortions in comparison with central schemes [6]. Among the various upwind discretizations available, Flux-Vector Splitting (FVS) and Flux-Difference Splitting (FDS) gained a large popularity in the last years. The idea of the first class is to decompose the convective flux into two parts according to the flow characteristics, as VanLeer [6]. Although their accurate representation of strong shocks and the excellent accuracy on supersonic flows, the resolution of boundary layers is quite poor due to their high numerical dissipation. Any improved version was proposed by Liou et al. [8] developing the so called Advective Upstream Splitting Method (AUSM) which is based on the decomposition of the inviscid flux into convective and pressure part. Even if the crisp resolution of strong shock and the accurate results of boundary layers, many authors have shown that spurious oscillations may be still present [8]. On other hand, the idea of the FDS is based on the solution of a locally 1D Euler equation at each face starting from a Riemann problem. In order to reduce the numerical effort required by the exact solution, approximated Riemann solvers were preferred, starting from the idea to use only partial information of the Jacobian structure. One of the most studied formulations is the modification of HLL to restore the missing contact surface, namely Harten-Lax-van Leer Contact (HLLC) in which the flux at the interface is computed by using three characteristic waves. This scheme is able to resolve both isolated shock and contact waves exactly under a suitable choice of contact and acoustic wave speeds. Unfortunately, this scheme suffers by numerical shock anomalies and postshock oscillations, especially under strong shocks [9]. Indeed, Godunov-type schemes fail to be accurate at low Mach numbers and they become unstable even for weakly compressible simulations. To overcome the previous issues, hybrid schemes, i.e. combinations of previous numerical classes, were investigated in order to provide a unique numerical scheme able to simulate almost every condition. This was motivated by Pandolfi et al. [9] which shown that such instabilities might occur because of mesh geometry, Mach number value and shock-capturing scheme. Hybrid schemes would then merge the different advantages of multiple discretization methods to possibly provide a unique robust and accurate numerical discretization. Particular attention was also given to the extension to quite low-Mach numbers in order to get a wide-range of applicability. Moreover, because to the excellent abilities already established of the HLLC scheme, the new version must include the HLLC scheme as reference scheme.

Hybrid scheme: HLLC-AUSM
The main idea was to make use of the excellent ability of the AUSM+Up-all speed scheme [8] on the resolution of boundary layer and low-Mach computation as well as its construction simplicity. To overcome the general issues of such scheme, HLLC scheme was employed. The idea is then to combine the HLLC scheme with the AUSM scheme to obtain a flux scheme capable to handle flows from low to quite high Mach number for steady and unsteady investigations. The inviscid flux ! " is split into two contributions, the convective term ! # and the pressure term ! $ : where the mass flux &̇ is discretized by using the classical HLLC formulation [10]. On other hand, the pressure term is discretized by the original AUSM+Up-all speed scheme. Finally, Schmidt et al. [11] proposed a similar hybridization using slightly different assumptions and testing it only on incompressible and subsonic applications.

Solution reconstruction and turbulence modelling
The previous scheme requires the information of the flow variables at the interface of the control volume requiring for a reconstruction method. The order of the spatial discretization is then determined by the accuracy of such reconstruction. To achieve at least a second order of accuracy in space in smooth regions, Monotone Upwind Schemes for Scalar Conservation Laws (MUSCL) [8] was employed in which the values stored at the centroid are used to estimate the values at the interfaces. More in details, the current implementation makes use of a piecewise linear reconstruction, where indicating left side (L) and right side (R) of the stored variable U, it claims: ) * = ) " + + " (∇) " ⋅ / * ) (2) ) 3 = ) 4 + + 4 5∇) 4

⋅ / 3 6 (3)
where + denotes the limiter function and 8 the distance from the centroid of the control volume. However, Blazek [6] shown that second-order upwind spatial discretization generates spurious oscillations in the regions where large gradients are present. The aim of the limiter function is then to restore the monotonicity of the solution in such region, in order to develop a monotonicity preserving scheme. The idea is to reduce the slope reconstruction to constrain the solution variations. In this implementation, two different slope limiters are employed, namely Barth-Jespersen [6] and Venkatakrishnan [6] following their classical definitions. Finally, Shock wave boundary layer interaction as well as separations and vortexes are common features which may occur into high-speed airbreathing engines. In order to reduce the computational cost, the turbulence terms are added implicitly as a source term, where for this first stage only URANS approach was considered because of its affordable computational cost.

Validation and Verification
The numerical implementation previously described was validated on multiple test cases, ranging from incompressible to supersonic, in order to establish its general accuracy. For showing the versatility, both structured and unstructured grids were employed to test the implementation as well as both Euler and Navier-Stokes models

Bidimensional cylinder
Steady-state inviscid simulation was performed around a cylinder of radius equal to 0.5m employing fully triangular unstructured mesh. Free stream inflow Mach number of 0.1 was imposed to test the splitting method in low-Mach simulation. The flow field, displayed in Fig. 1, shows the typical symmetric velocity contours distribution. Moreover, the pressure coefficient of the upper side was compared with the analytical one depicting an excellent agreement.

Laminar boundary on a flat plate
The accuracy of the numerical implementation was tested simulating the boundary layer over an adiabatic flat plate of length L. Structured mesh was employed to perform such analysis, using a y + =1.

Bidimensional NACA 0012 airfoil
To demonstrate the ability of simulating flows where both subsonic and supersonic regions coexists (transonic regime), the NACA 0012 airfoil was simulated at zero angle of attack imposing a free stream Mach equal to 0.8 employing fully unstructured mesh. The qualitative results in terms of pressure contour and streamline are reported in Fig. 3. The supersonic pocket is clearly visible on the suction side where it ends up with a normal shock. As well in Fig. 3 is presented also the trend of the pressure coefficient, on the suction side, compared with the experimental data of Farrokhnal et al. [12]. The trend is well reproduced comparing it with the experimental data, where the slight difference is due to the inviscid flow assumption.

Supersonic backward facing step
The backward facing step is a classical validation test case in which expansion fan, recirculation region and weak oblique reattachment shock are present at the same time. The inflow Mach number of 2.5 was imposed at the inlet and k-w SST was employed to solve the turbulent effects. The usual recirculation zone after the step is displayed in Fig. 4 with streamlines details as well. According to the measurements performed by Smith et al. [13], the recirculation zone should take place around 2.4h leading to a percentage error less than 10%. Moreover, the pressure after the step was compared with Smith's experimental data (Fig. 4). The hybrid HLLC/AUSM scheme shows a perfect agreement with the experimental position leading to a really precise reattachment location with respect to other solvers which fail such feature. Must be noted that a proper tuning of the boundary conditions of the turbulence model was needed to increase the quality of the validation.

Fuel-off 2D scramjet engine
The overall implementation was finally verified on a test case available in literature, typical for propulsion applications. The 2D scramjet geometry was investigated both numerically and experimentally by Lorrain et al. [7]. The inlet consists of two ramps of different slopes, to promote the compression, leading to a constant area combustor which ends with a nozzle of constant expansion ratio. The combustion chamber part was simulated in fuel-off condition. This test aims to demonstrate the ability of the solver to handle different flow features at the same time, as strong expansions, shock reflections, shock wave-boundary layer interactions and expansion fans. A freestream Mach number of 7.5 was imposed at the inlet and a structured mesh was built in order to avoid the use of wall functions. Experimental investigations were conducted by Lorrain et. al [7]. A qualitative comparison about the inlet area is presented in Fig. 6 in terms of density gradient (Schlieren). The shock wave nature of the flow is well represented as well as the locations. More fundamental, the reflections are still in good agreement where however, some narrow small instabilities are present just after the corner. Must be mentioned that, in this case, the dissipative nature of the turbulence and the high-aspect ratio of the grid could drastically change the solution.

Conclusion
An advanced open-source compressible solver was developed using openFoam toolbox for providing high-fidelity simulations of combined propulsion systems for hypersonic  (Smith, 1967) aircrafts. Stability, accuracy and low computational cost were taken as the reference requirements. On this view, customizable explicit time integration scheme suitable for both unsteady and steady investigations was implemented. An innovative hybrid scheme for inviscid flux discretizations was carried out following the idea of merging AUSM and HLLC implementations. MUSCL reconstruction was taken into account to increase the spatial accuracy by using two different multi-dimensional flux limiters. At the end, standard viscous flux discretization as well as turbulence models were linked to solve turbulent compressible Navier-Stokes equations. The whole implementation was then validated in various test cases, ranging from incompressible to hypersonic supersonic regimes, at different flow conditions. The new inviscid flux discretizations shown a good agreement in the case of incompressible where, on supersonic test cases, the hybrid scheme shown a better ability to capture expansion phenomena, shock wave-boundary layer interaction and contact discontinuities.