%~Mouliné par MaN_auto v.0.37.1 (6604bb9b) 2025-12-10 13:43:19
\documentclass[CRMECA,Unicode,biblatex]{cedram}

\addbibresource{CRMECA_Lugrin_20250928.bib}

\usepackage{subcaption}
\usepackage{siunitx}

\newcommand{\dis}{\mathrm{dis}}
\newcommand{\tot}{\mathrm{tot}}
\newcommand{\wall}{\mathrm{wall}}

\DeclareMathOperator{\Rey}{Re}

\makeatletter
\renewcommand\p@subfigure{}
\makeatother

%%%---------------------------------------------------------------------------

%%% DELIMITERS

\usepackage{mathtools}
\DeclarePairedDelimiter{\parens}{\lparen}{\rparen}
\DeclarePairedDelimiter{\abs}{\lvert}{\rvert}
\DeclarePairedDelimiter{\bracks}{[}{]}

%%%---------------------------------------------------------------------------

\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\title{Towards an automated toolset for mesh adapted Navier--Stokes simulation of hypersonic vehicles}
\alttitle{Vers un ensemble d'outils automatisés pour la simulation Navier--Stokes avec adaptation de maillage des véhicules hypersoniques}

\author{\firstname{Mathieu} \lastname{Lugrin}\IsCorresp}
\address{DAAA, ONERA, Institut Polytechnique de Paris, 92190 Meudon, France}
\email{mathieu.lugrin@onera.fr}

\author{\firstname{Baptiste} \lastname{Isnard}\IsCorresp}
\address[1]{DAAA, ONERA, Institut Polytechnique de Paris, 92190 Meudon, France}

\author{\firstname{Clément} \lastname{Benazet}}
\address{DAAA, ONERA, Institut Polytechnique de Paris, 92320 Châtillon, France}

\author{\firstname{Bruno} \lastname{Maugars}}
\address[3]{DAAA, ONERA, Institut Polytechnique de Paris, 92320 Châtillon, France}

\author{\firstname{Cédric} \lastname{Content}}
\address[3]{DAAA, ONERA, Institut Polytechnique de Paris, 92320 Châtillon, France}

\thanks{This project was provided with computer and storage resources by GENCI thanks to the grant 2023-A0142A14106 on the SKL partition of the supercomputer Irene at TGCC. The development of the SoNICS CFD software used for this study was partly funded by the DGAC SONICE project}
\CDRGrant[GENCI]{2023-A0142A14106}

\keywords{\kwd{Hypersonic} \kwd{RANS} \kwd{mesh adaptation}}
\altkeywords{\kwd{Hypersonique} \kwd{RANS} \kwd{adaptation de maillage}}

\editornote{Article submitted by invitation}
\alteditornote{Article soumis sur invitation}

\begin{abstract}
A fully automated CAD-to-post toolset for Navier--Stokes and Reynolds Averaged Navier--Stokes (RANS) simulations of high speed vehicles flow is presented. The toolset combines the vertex-centered SoNICS solver with anisotropic mesh adaptation based on the REFINE toolbox, eliminating manual meshing and enabling efficient, parallelized simulations from CAD input to converged solutions. The tool is validated on open cases representative of the complexity of the flow encountered around high-speed vehicles including reentry and airbreathing cruise vehicles. This includes complex three-dimensional forebody laminar simulations and axisymmetric triconic cases. For the internal aerodynamics part, given the lack of proper validation references, a direct numerical simulation of the internal conduit of a generic dual mode ramjet air intake is conducted and compared to legacy RANS simulation and mesh-adapted RANS simulation. A demonstration on a full scramjet-powered cruise vehicle then highlights the toolset’s scalability and adaptability to industrial applications. The automated workflow reduces setup time and human error, making it a valuable asset for hypersonic vehicle design and optimization.
\end{abstract}

\begin{altabstract}
Un ensemble d'outils entièrement automatisé pour les simulations Navier--Stokes et Reynolds Averaged Navier--Stokes (RANS) des écoulements autour des véhicules à grande vitesse est présenté. Cet ensemble d'outils combine le solveur SoNICS avec l'adaptation de maillage anisotrope basée sur la boîte à outils REFINE, éliminant ainsi le maillage manuel et permettant des simulations parallélisées efficaces, depuis l'entrée CAO jusqu'aux solutions convergées. L'outil est validé sur des cas ouverts représentatifs de la complexité de l'écoulement rencontré autour des véhicules à grande vitesse, y compris les véhicules de rentrée et les véhicules de croisière à propulsion aérobie. Cela inclut des simulations laminaires tridimensionnelles complexes de l'avant du fuselage et des cas triconiques axisymétriques. Pour la partie aérodynamique interne, compte tenu du manque de références de validation appropriées, une simulation numérique directe du conduit interne d'une prise d'air générique pour statoréacteur est réalisée et comparée à une simulation RANS classique et à une simulation RANS avec adaptation de maillage. Une démonstration sur un véhicule de croisière propulsé par statoréacteur souligne ensuite l'évolutivité et l'adaptabilité de l'ensemble d'outils aux applications industrielles. Le flux de travail automatisé réduit le temps de configuration et les erreurs humaines, ce qui en fait un atout précieux pour la conception et l'optimisation des véhicules hypersoniques.
\end{altabstract}

\begin{document}
%\input{CR-pagedemetas}
%\end{document}
\maketitle

\vspace{3\baselineskip}

\section{Introduction}

The design of high-speed vehicles heavily relies on Computational Fluid Dynamics (CFD), primarily through steady Navier--Stokes (NS) or Reynolds-Averaged Navier--Stokes (RANS) simulations. These simulations provide access to key performance indicators such as integrated data (e.g., lift to drag ratio, thrust, aerodynamic stability, etc.) or local quantities (heat fluxes, pressure load, etc.) that are then used to evaluate and refine the designs. Even with recent advances in terms of tools and solvers, setting up a CFD simulation from a CAD drawing of a system is still a time-intensive task that must be performed by a skilled engineer. While the use of unstructured meshes has significantly reduced the human time required to build medium to high quality meshes, meshing remains one of the main blocking points towards the use of CFD in automated toolchains. As the aerospace industry continues to push for more efficient design, the necessity for automated CAD to post-process workflows becomes increasingly evident. The manual CAD to mesh process, which is prone to human error and time-consuming iterations, obstructs the overall productivity of engineering teams. An integrated, automated meshing solution would significantly streamline the CFD pipeline, enabling faster design evaluation and optimization.

A possible solution for automatic CAD to post from the literature is the use of Immersed Boundary Method (IBM) \cite{mittal2005immersed}. The main advantage of this technique is that it does not require a body fitted mesh, and thus can use automatically generated Cartesian or octree meshes coupled with automatic masking of the area inside the solid. This method has been widely used for lattice-Boltzmann and steady or unsteady Navier--Stokes simulations, such as wall-modelled LES~\cite{verzicco2023immersed}. However, it usually suffers from multiple drawbacks. First, due to the absence of body fitted meshes, IBM simulations often rely on wall models, which have the common issue of lacking grid convergence~\cite{terracol2025dynamic} and exhibiting spurious behavior in separated regions~\cite{shih2003application}. This last drawback may be problematic for high-speed vehicle evaluation, particularly for cruise vehicles, where the air intake uses multiple shock-boundary layer interactions to slow down the flow and often reaches a fully separated state. Getting precise laminar baseflow using wall-model can also be challenging. Then, the use of body forces inherently leads to a non-conservative solution, which can again be a major problem for internal aerodynamics, as the mass imbalance may cause a wrong evaluation of the engine performance.

Another very promising solution to alleviate the manual meshing problems without the drawbacks of IBM is anisotropic mesh adaptation. It provides good convergence of the solution even with no physical intuition nor knowledge in the starting mesh, and do not suffer from conservativity or spurious behavior in separated regions, given that it remains resolved and body fitted. Anisotropic mesh adaptation is known for its high efficacy in the automatic generation of optimal computational meshes, first within the context of inviscid flow simulations~\cite{Loseille10,ALAUZET2010561,loseille2007,VENDITTI200240} and more recently for RANS simulations~\cite{frazza2021,MORISCO202557}. At its heart, mesh adaptation~\cite{Castro1997,Dervieux2003,FREY20055068,GRUAU20054951,Li2003} involves an iterative process: utilising data from a preceding simulation to automatically dictate the necessary mesh adjustments required to enhance either the global solution accuracy or a particular quantity of interest. This is achieved via solution-based sensors that analyse the current results and subsequently define modifications to the mesh. The goal of these prescribed modifications is to balance and optimise the element distribution ---~and, consequently, the computational resources~--- throughout the domain~\cite{Fidkowski2011a,Roge2008,Hecht2008}. Regions characterised by significant gradients in the solution, or those deemed important for the computation of a specific output, are selectively and automatically refined to increase the local accuracy of the discretization. The validity and optimality of this methodology are underpinned by a formal mathematical analysis of the associated numerical error~\cite{LoseilleHDR}.

However, most toolchains relying on mesh adaptation are still not automated CAD-to-post workflows, as the tool often requires manually created surface and initial mesh of the model as well as the tagging of the boundary conditions. They are also often relatively inefficient, even if the mesh adaptation by itself allows for more efficient simulation due to the lower amount of points needed to converge a solution, as they chain tools through Input/Output (I/O, here used to refer to the writing and reading of the solution on a physical disk space) and different executables. Most of the tools are also not built to work in a parallel or distributed manner, which strongly limits performance and, worse, the maximal number of points that can be computed, as one may reach the memory limit of the computing nodes.

\goodbreak

In this context, this paper presents a fully automated, parallel distributed CAD-to-post toolset for NS and RANS simulation based on anisotropic mesh adaptation, with a particular focus on the demonstration of its ability to capture the flow around high-speed vehicles.

While the tool has already been extensively validated during development using non-regression tests from the elsA solver, as well as open literature databases such as the NASA Turbulence Modeling Resource, conducting simulations of high-speed vehicles introduces new complexities that are specific to hypersonics and may not be covered by the usual CFD code validation. The tool must be able to handle highly compressible flow, including strong shocks, shock-shock, and shock-boundary layer interactions. This may also be challenging in terms of mesh adaptation, given that metric based adaptation will tend to refine regions with the strongest gradients (i.e.\ shocks) while overlooking other physically important regions of the flow, such as boundary and mixing layers.

To assess the ability of the tool to capture the flow around high-speed vehicles, it is tested on open cases representative of the complexity of the flow encountered in real-life applications, namely: a canonical reentry vehicle, a complex 3D forebody, and a supersonic internal air intake. Finally, a demonstration on a mock-up full cruise vehicle is conducted.

The article is organised as follows. First, the Navier--Stokes solver, the re-mesher, and the way they are linked together are presented. The manner in which the initial geometry, first mesh generation, and boundary condition tagging are handled is also highlighted, given that they are of great importance for the automation process. A new treatment of the metric used for adaptation and CAD handling is also presented. The tool is then validated on different open cases from the literature, representative of the complexity of real life applications.

The first case of interest is a canonical reentry vehicle that has been widely studied both numerically and experimentally~\cite{caillaudAIAAJ,benitez2021separation}. It consists of an axisymmetric cone-cylinder-flare geometry, which offers numerous advantages. First, it allows for the comparison of the results and performances of the toolset on a bi-dimensional case against well established solvers from the literature. Even if geometrically simple, this case includes most of the complexity of the flows encountered around reentry vehicles.

The second and third cases of interest relate to airbreathing cruise vehicles. First, an external aerodynamic case is needed to validate that the tool can handle hypersonic forebodies, potentially with complex flow features such as curved shocks, the linked vortices, and high-Reynolds number laminar boundary layer. One main point of interest here will be the ability of the tool to capture accurate wall heat fluxes, which is both a challenging quantity to obtain numerically and one that is of paramount importance for the design of vehicles, as it will constrain the thermal protection system. Here, the BOLT case is chosen as a reference given that it provides geometrical and physical complexity close to that of real vehicles while being open and widely studied in the literature, with experimental (both ground~\cite{sombaert2025ground,chynoweth2025boundary} and flight~\cite{wheaton2018boundary,wirth2024measurements}) and numerical results~\cite{vu2024postflight}. Specifically, a very challenging numerical reproduction of the highest laminar flight Reynolds number of the  BOLT~II campaign is conducted and compared to numerical results from~\cite{vu2024postflight}.

The last validation case of interest relates to internal aerodynamics linked with ramjet propulsion. To get a properly functioning ramjet, one needs to slow down and compress the incoming supersonic flow to reach a subsonic state inside the engine. This can be done using both external and internal compression apparatus. Using a series of shock-waves, the flow is slowed down to feed the combustion chamber. To properly design the engine, one needs to have precise data about the flow that will enter the combustion chamber. Being able to easily and accurately compute such flows is thus mandatory to correctly predict vehicle trajectories and validate designs. For this case, the objective is to demonstrate that the tool can reliably reproduce the key internal flow features, such as corner vortices, shock boundary layer interaction (SBLI) and shock-train evolution. Even if some canonical cases of shock trains are available in the literature~\cite{fievet2016numerical,carroll1988numerical}, they are too far from applicative needs to be used in the scope of this paper. As such, a new validation case was built from scratch to reproduce all the flow features of interest. To establish a reference, a DNS of the benchmark configuration is performed, providing high-fidelity data against which results from the tool is compared. Given that the tool relies on turbulence modelling to capture such flow, a ``legacy'' RANS simulation (conducted with the elsA~\cite{elsA} well established solver) was conducted with the same models in order to discriminate discrepancies between the DNS and the RANS simulations caused by the turbulence modelling from those that may be caused by an ill behaviour of the toolset. Finally, once the ability of the tool to capture both complex internal and external flow has been proven, a demonstration on a full mock-up vehicle (with a modelled scramjet engine) is presented as a demonstration of the possible use of the tool.

\section{Toolset description}

The toolchain (presented in Figure~\ref{fig:toolchain}) is built based on two main tools, namely the vertex-centered Navier--Stokes solver SoNICS from ONERA--Safran~\cite{sonics} briefly described hereafter and an anisotropic mesh adaptation toolbox adapted from the open-source toolbox REFINE developed by NASA~\cite{REFINE_verification}. Every step that will be described in the following is linked through memory APIs. This means that once the CAD file have be read by process, no I/O has to be conducted (except for post-treatment needs). Also, all following steps are executed in a parallel context, but parallel management is completely transparent to the user.

\subsection{CFD Solver}

SoNICS is a CFD software that solves Euler, Navier Stokes or RANS equations using either cell-centered or vertex-centered second-order finite volume schemes. It can perform multi-species reactive flows simulations. It has been designed with four major requirements: it must be easy to configure and must be able to take advantage of powerful hardware (CPU and GPU), which are all complex, distributed and heterogeneous~\cite{DOEEexascaleBES}; it should handle a large variability in terms of Riemann solvers~\cite{roe_jcp_1981,CHEN2018126,Renac2021}, flux limiters~\cite{VENKATAKRISHNAN1995120,VANALBADA1982,PIPERNO1998183}, fluid models (equations of state, viscosity laws, mixing laws...), turbulence modelling (RANS, LES, DES); it should be fully differentiable in order to allow modern optimization algorithms~\cite{kenway2019}, input-output analysis~\cite{Trefethen1993,Alizard09,McKEON_SHARMA_2010,Sipp2013,Rigas_Sipp_Colonius_2021,bugeat20193d}, multi-disciplinary optimization~\cite{Martins2013,Kenway2014}, data assimilation~\cite{10.1016/j.jcp.2015.11.012,Duraisamy2019,LI2017111,PhysRevFluids.5.094603}, end-to-end machine learning~\cite{NEURIPS2019_01386bd6}; and finally, it should provide APIs to help couple the tool with other solvers. A full description of the architecture in use for the SoNICS solver can be found in the companion paper by Lienhardt and Michel~\cite{lienhardt}.

As mentioned above, the vertex-centered solver of SoNICS is used. The spatial discretization is based on the one proposed by A.~Dervieux and C.~Debiez~\cite{DEBIEZ2000} on simplex mesh (triangle/tetrahedron). The convective terms are solved by the finite volume method on the dual mesh composed of median cells. Fluxes are determined through the use of the AUPM approximate Riemann solver proposed by Chen et~al.~\cite{CHEN2018126}. Second order space accuracy is achieved through MUSCL type reconstruction method~\cite{VANLEER1997}. The MUSCL procedure is applied as proposed in A.~Dervieux and C.~Debiez~\cite{DEBIEZ2000} or Alauzet and Frazza~\cite{ALAUZET2021110340} to recover the so-called V4-scheme. A shock capturing similar to the one of~\cite{sciacovelli} is also implemented. Viscous fluxes are computed with the conventional expression described as Cell Based Viscous scheme in Liu et~al.~\cite{Liuetal2023}.

Boundary conditions are weakly imposed by evaluating the flux at each facet at the boundaries of the domain. As in Alauzet and Frazza~\cite{ALAUZET2021110340}, each Dirichlet boundary (such as the velocities on a no-slip wall) are strongly enforced at each iteration, i.e.\ the DOF concerned by the boundary are removed from the whole CFD problem.

In order to reach the steady state, we consider an implicit backward Euler scheme with local time stepping, approximately solved by an efficient LU-SGS solver similar to the one proposed in~\cite{SONG2022}. For RANS simulation, the Spalart--Allmaras (SA) \cite{SPALARTAL} model with the compressibility correction of~\cite{DeckSA2002} and the Quadratic Constitutive Relation (QCR) \cite{spalart2000strategies} are used.

\subsection{Anisotropic mesh adaptation on full tetrahedron meshes workflow}

\begin{figure}[h!]
\includegraphics[width=\textwidth]{workflow.png}
\caption{Diagram of the toolchain used for the automated computations.}\label{fig:toolchain}
\end{figure}

The automated toolchain for mesh adaptation is presented in Figure~\ref{fig:toolchain}. The input of the process is a CAD description of a closed shell representing the fluid volume to simulate. The first step is then to prove the existence of the mesh. For this purpose, REFINE provides an EGADS~\cite{EGADS} wrapper to generate a surface mesh of the CAD shell. One can note that for the applications presented here, the CAD model is provided as a STEP file, which allows interaction with any CAD modeling software; however, a direct link with ESP~\cite{ESP} could also be utilized to take advantage of its differentiation. From the STEP input, EGADS tessellates the CAD edges and surfaces to build a closed surface mesh. Tesselation provides meshes whose elements have poor quality for a volume mesher, especially given that each CAD patch is meshed independently. As such, the tesselation process results in an overall low-quality surface mesh. In order to improve its quality before providing it to the tetrahedral mesher, REFINE performs surface mesh adaptation driven by a curvature-based metric that is evaluated thanks to the CAD parametrisation at each mesh node. Several passes can be performed with various parameters (minimal and maximal cell size, mesh complexity, size gradation) in order to obtain a decent surface mesh before volume mesh generation. Once done, a tetrahedral mesh generator (here tetgen~\cite{hang2015tetgen} or netgen~\cite{schoberl1997netgen}) is used to mesh the interior of the fluid volume. This first mesh generation is illustrated in Figure~\ref{fig:firstmesh}.

\begin{figure}[h!]
\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{image (3).png}
\caption{}\label{fig:firstmesh_a}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{image (1).png}
\caption{}\label{fig:firstmesh_b}
\end{subfigure}

\vskip .5 \baselineskip

\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{image (2).png}
\caption{}\label{fig:firstmesh_c}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{image (4).png}
\caption{}\label{fig:firstmesh_d}
\end{subfigure}
\caption{Automated first mesh generation process illustrated by: \eqref{fig:firstmesh_a}~CAD model; \eqref{fig:firstmesh_b}~the EGADS tesselation of the CAD model; \eqref{fig:firstmesh_c}~the surface mesh after remeshing by REFINE; \eqref{fig:firstmesh_d}~volume mesh created by tetgen.}\label{fig:firstmesh}
\end{figure}

The mesh structure ---~including coordinates, connectivity, and CAD information~--- is converted into a Python/CGNS structure~\cite{rumsey2012recentCGNS} (CGNS tree), which serves as the entry point for the SoNICS solver. The CGNS tree data is natively distributed over MPI processes following Maia's parallel rules for CGNS trees (see~\cite{maia,maia_doc}). It will remain distributed for the rest of the workflow.

In industrial configurations, the number of ridge and surface patches can be exceedingly large, which implies some extra work for the user to retrieve link between CAD patches and the CFD boundary conditions (BCs). In order to relieve the user of this step, if CAD patches are tagged with names using the advanced face option of STEP format, the toolchain preserves these names during the mesh generation process, so boundaries can easily be retrieved for the solver boundary condition tagging. Patches from the same families can also be concatenated together in a single CGNS node to ease I/O and tree handling.

The Python/CGNS tree is then prepared by adding a family of boundary conditions (BCs) (automatically, thanks to the preservation of patch names), the reference state, and the solver options. Steady compressible Navier--Stokes or RANS simulations are then conducted using the vertex-centered solver of SoNICS described here-before. Once a solution is found, the computation of the mesh adaptation metric is performed directly in the CFD solver and is fed to the remesher. In this work, the metric is computed as the absolute value of the Hessian of the Mach number (albeit the user can decide to use any field). A pre-processing is applied by REFINE to compute the error estimate for the desired $p$-norm and complexity. Then an anisotropic size gradation algorithm is used to smooth the metric. REFINE then performs a parallel remeshing of the surface and volume meshes. The geometric description of the fluid domain is preserved with high fidelity using EGADS's CAD projection when inserting or moving boundary mesh vertices. The current flow solution is finally interpolated on the new adapted mesh using the open-source ParaDiGM library~\cite{paradigm} (see also~\cite{andrieu2025}), which is available through the Python/CGNS API of Maia. The interpolated solution is used as an initial state and the simulation is run again. The process is looped, potentially with an increase in complexity depending on the steps, until a proper solution is converged. Convergence indicators are case-dependent and are to be included in the Python script directly by the person skilled in the art. The same goes for complexity jumps, gradation, and most of the CFD solver parameters which are also application-dependent.

The whole procedure that is described above is performed by a single Python script that takes as inputs the STEP file and the run conditions. All pre/co/post-processing, partitioning and distribution are performed using the open-source Maia library.

\subsection{Convergence}

Given that the first generated mesh is very rough and does not include any physical intuition, the first NS or RANS computation at the beginning of the loop can be numerically challenging. The usual way to start the mesh convergence process is to run the solver computation over a few pseudo-time iterations to avoid solver crashes. The number of iterations can then be slowly increased over remeshing passes to gradually converge the proper solution. A trick allows to accelerate this solution convergence process: the Mach Hessian used as a metric for feature-based remeshing is a good criterion for mesh convergence of Navier--Stokes equations resolution~\cite{loseille2011continuouspartII}. However, if used as is, it can lead to slow mesh convergence because it is not very efficient at quickly capturing the boundary layer. Especially for flows including shocks, in the first steps where meshes have poor discretization, the chosen metric highlights the very strong gradients in compression regions compared to viscous regions. Frazza~\cite{frazza2021} proposed to boost the Hessian near viscous walls in order to add more points into the boundary layer, which will help the process converge quickly to a steady flow solution. For this purpose, metric diagonal is boosted by a factor in the wall normal direction.

\section{Reentry vehicles aerodynamics}

As stated in the introduction, the first type of high-speed vehicles of interest are the reentry vehicles. To assess the ability of the solver to treat those cases, the axisymmetric cone-cylinder-flare CCF12 geometry, which was designed as a canonical and open case to study the flow features found in atmospheric reentry vehicles is chosen as a benchmark. It includes a simple developing hypersonic laminar boundary layer followed by more complex features such as boundary layer interaction with expansion fans, and most importantly, separation induced by the shock-boundary layer interaction. The numerical results of Caillaud et~al.~\cite{caillaudAIAAJ} (see Table~\ref{tab:flow-conditions1}) are used as a reference given that they provide high-order solutions cross-validated between two well established solvers from ONERA~\cite{poulain2023broadcast} and NASA~\cite{baurle2022vulcan}. For this case, the axisymmetric solver of SoNICS is used alongside the 2D remesher of REFINE, in this case the mesh is only composed of triangles.

\begin{table}[b]
\caption{Flow conditions from~\cite{caillaudAIAAJ} used in this study.}\label{tab:flow-conditions1}
\resizebox{\textwidth}{!}{
\begin{tabular}{cccccccccc}
\hline
\\[-2.7ex]
\hline
\\[-2.5ex]
$\Rey_m \times10^6$ (\unit{\per\meter}) & $R_n$ (mm) & $u_\infty$ (m/s) & $\rho_\infty$ (kg/m$^3$) & $T_\infty$ (\unit{K}) & $T_0$ (\unit{K}) & $T_{\wall}$ (\unit{K}) & $P_0$ (kPa) \\ 
\hline
\\[-2.7ex]
\hline
\\[-2.5ex]
2.229 & 0.1 & 907.92 & 0.00920 & 56.977 & 467.21 & 300 & 237.50 \\
4.614 & 0.1 & 990.41 & 0.02128 & 67.800 & 555.96 & 300 & 653.90 \\
\hline
\end{tabular}}
\end{table}

In that specific 2D axisymmetric case, the CAD file is generated from the analytical description of the body using GMSH and exported as a STEP file. An example of an adapted mesh on that case, alongside a Mach number contour, is presented in Figure~\ref{fig:mesh_ccf}.

\begin{figure}[h!]
\includegraphics[width=0.7\textwidth]{mesh_ccf.png}
\caption{Adapted mesh and Mach number contour for the low Reynolds number case of~\cite{caillaudAIAAJ}.}\label{fig:mesh_ccf}
\end{figure}

Taking the low Reynolds case as an example, the toolset is run for 121 seconds on 4 cores of an Intel Cascade Lake CPU, allowing us to conduct 6 mesh adaptations and steady computations until complete convergence is reached, the evolution of the wall pressure distribution for each step is presented in Figure~\ref{fig:conv}.

\begin{figure}[h!]
\includegraphics[width=0.48\textwidth]{axi_conv.png}
\caption{Evolution of the wall pressure distribution on the CCF12 at low Reynolds number at each iteration showing the convergence path.}\label{fig:conv}
\end{figure}

The wall pressure distribution for two different Reynolds numbers is presented in Figure~\ref{fig:valid_ccf12} alongside results from two high-order solvers. Given the sensitivity of the strongly separated shock boundary layer interaction to both the resolution of the boundary layer gradients and the associated separation and reattachment shock, obtaining the correct pressure signature guarantees that these features are accurately captured in the solution. This result validates the good behaviour of the toolset for hypersonic axisymmetric configurations.

\begin{figure}[h!]
\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{axi_22.png}
\caption{Low Reynolds case}%\label{subfig:evp_convergence}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{axi_46.png}
\caption{Intermediate Reynolds number case}%\label{subfig:res_convergence}
\end{subfigure}
\caption{Wall pressure distribution for the two different Reynolds number alongside reference data from~\cite{caillaudAIAAJ}.}\label{fig:valid_ccf12}
\end{figure}

\section{Cruise vehicles aerodynamics}

The second type of vehicle of interest is air-breathing cruise vehicles, which, in addition to external aerodynamics, involve complex internal flow. This validation is split into two parts: first, a 3D forebody case is studied; then, a canonical internal aerodynamic reference case is built and used as a reference to validate the toolset.

\vskip 1.5\baselineskip

\subsection{Forebody}

\vspace{\baselineskip}

\begin{table}[b]
\caption{Flow conditions from~\cite{vu2024postflight} used for the BOLT forebody reference, data coming from the BOLT-II flight.}\label{tab:flow-conditions}
\resizebox{\textwidth}{!}{
\begin{tabular}{lccccccccc}
\hline
\\[-2.7ex]
\hline
\\[-2.5ex]
Time (s) & $\Rey_\infty \times10^6$ (\unit{\per\meter}) & Altitude (m) & $u_\infty$ (m/s) & $\rho_\infty$ (\unit{kg/m^3}) & $T_\infty$ (\unit{K}) & Mach & $T_{\wall}$ (\unit{K}) \\
\hline
\\[-2.7ex]
\hline
\\[-2.5ex]
401.9 & 6.0914& \num{23908} & 1768.5 & 0.048302 & 213.15 & 6.0424 & 390 \\
\hline
\end{tabular}}
\end{table}

The first subcase of interest for cruise vehicles concerns forebodies, which are often complex 3D geometries on which the main points of interest will be drag, lift and heat-flux prediction. The BOLT case is chosen as it provides plenty of open data (including flight) and its 3D geometry is representative of the complexity of real designs. It also displays an intricate baseflow structure with multiple vortices, which is challenging to capture and is thus a good benchmark. The chosen condition for the computation is a challenging high Reynolds case with both CFD and sparse flight data coming from~\cite{vu2024postflight}. This case is at a flight altitude where laminar-turbulent transition starts to appear, which means that it is the highest Reynolds number for which the flow remains laminar; the full flow conditions are presented in Table~\ref{tab:flow-conditions}. Given the high Mach number and the relatively low altitude, this case is fully representative of a scale one hypersonic cruise vehicle forebody in flight. This case would also serve future validation of transition prediction tools integrated to the solver.

The case is run on a quarter of the full scale flight  BOLT~II model using symmetry planes. The computation is conducted on 2400 Intel Cascade Lake cores for a single job of 15 hours, allowing for 8 adaptation/computation steps and achieving a final mesh of 61~million nodes in order to reach proper convergence of the wall heat fluxes. This can be compared to the 110M~cells used in the study of~\cite{vu2024postflight}, keeping in mind that the solution on the structured grid is computed using a sixth order central finite difference scheme, while the mesh adapted solution uses a second order method. This gives insight on the ability of mesh adaptation to strongly reduce the required number of DOF to achieve the same convergence. This could be even more drastic when using goal-orientated adaptation (on the parietal heat-flux) \cite{Puigt2025MeshAdaptation}, given that the convergence of wall heat-fluxes requires a lot of degrees of freedom, especially with meshes composed only of tetrahedron and using metric based adaptation.

An illustration of the geometry, alongside the wall heat-flux distribution and a Mach number field at the outlet of the domain, is presented in Figure~\ref{fig:overviewbolt}. It highlights the complex nature of the flowfield with multiple vortical structures appearing, especially in the center region. Quantitative comparisons of heat-fluxes at two different longitudinal positions with the CFD solution of~\cite{vu2024postflight} are presented in Figure~\ref{fig:boltslices}. Overall, excellent agreement is found between the two simulations. Results from the mesh adaptation display a relatively good smoothness given that they come from anisotropic tetrahedron only meshes. Some discrepancies are still visible in the peaks linked with local vortices, mostly in terms of amplitude of the fluctuation. Figure~\ref{fig:boltslices}\eqref{subfig:evp_convergence} presents the most upstream slice, from left to right, the relatively low heat-flux region is linked with the upwash of the centerline vortices (caused by the 3D nature of the geometry and thus of the bow shock which leads to the creation of vorticity), those vortices are clearly visible in Figure~\ref{fig:overviewbolt}. This region is then followed by a first peak linked with the downwash of those vortices and then a long plateau of relatively high heat flux linked with the thin boundary layer in the outboard region. In all those regions, the two computations are in almost perfect agreement.

The origin of these discrepancies is still unknown; they may be linked to differences in the mesh resolution; some vortices may be under-resolved in either simulation. Another possible origin could be the growth of vortices linked to numerical errors from the time integrators, as the baseflow is extremely unstable; small numerical errors could lead to the growth of structures that are not part of the baseflow. Future work on that case will require additional simulations from other meshes and/or solvers to assess the origin of these discrepancies. Again, a Newton method could be used once the Jacobian operator is available to compute a fixed point of the equations.

Overall, the tool appears to be fully capable of computing a very complex baseflow such as the one encountered on the BOLT forebody. The key point to highlight here is that the user time needed to set up such a computation is minimal (of the order of minutes) with no effort whatsoever put into meshing, which strongly contrasts with other studies on that case.

\begin{figure}[h!]
\includegraphics[width=1\textwidth]{visu_bolt.png}
\caption{BOLT geometry and simulation results (wall heat flux and Mach number field at the outlet of the simulation domain) for the descent case of~\cite{vu2024postflight}.}\label{fig:overviewbolt}
\end{figure}

\begin{figure}[h!]
\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{0.3579.png}
\caption{}\label{subfig:evp_convergence}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{0.796.png}
\caption{}\label{subfig:res_convergence}
\end{subfigure}
\caption{Spanwise heat-fluxes distribution at two longitudinal sections of the BOLT forebody alongside data from the US3D~\cite{candler2015development} computation of~\cite{vu2024postflight}.}\label{fig:boltslices}
\end{figure}

\subsection{Internal aerodynamics case}

As explained in the introduction, given the lack of open cases representative of applied vehicle technologies, a reference case is built to validate the toolset on hypersonic internal aerodynamics. The motivation for the geometry and flow conditions choices as well as the setup for the reference simulation is presented in this section, followed by comparison between reference, legacy RANS and automatic toolset results.

\subsubsection{Geometry and flow conditions}

The geometry under study is designed to be representative of the complexity of the internal flow of a ramjet powered vehicle. As such, it should include an internal compression region, followed by a throat, and then a slowly diverging section in which a shock train will appear. A back pressure at the outlet of the diverging section is necessary to simulate the engine. To facilitate the sharing and accurate reproduction by other researchers, the geometry must be canonical and analytically describable. This led to the geometry described in Appendix~\ref{app:geom}. The dimensions of the geometry, while arbitrary, are chosen such that reproductions in conventional supersonic wind tunnels are facilitated. The flow conditions are chosen to match a reasonable value for most supersonic wind tunnels to allow for further cross comparison.

The Reynolds number is chosen as high as possible, while still allowing for a highly resolved computation within the given computing allocation. The inlet Mach number is set to 3, the stagnation temperature is low (300~\unit{K}) to facilitate experimental reproduction, and the static pressure of the inlet is chosen to adjust the Reynolds number and thus the computational cost. The complete conditions, as well as the reference values used for the computations, are presented in Table~\ref{tab:table1}.

\begin{table}[hbt!]
\caption{Reference conditions for the computations.}\label{tab:table1}
\begin{tabular}{cccccccc}
\hline
\\[-2.7ex]
\hline
$M_\infty$ & $P_\infty$ & $T_i$ & $ \alpha$ & $\beta$ &$T_{\wall}$ & $\Rey_L$ & $L$\\
\hline
\\[-2.7ex]
\hline
\\[-2.5ex]
3 & 4000 Pa & 300 K & \ang{0} & \ang{0} & 300 K & $4.029\times10^6$ & 0.4 m\\
\hline
\end{tabular}
\end{table}

The boundary conditions are supersonic inlet and outlet, and isothermal viscous wall at 300~K. The turbulence generation technique is described in Appendix~\ref{app:turbu}. The chosen conditions lead to the approximate Reynolds number presented in Table~\ref{tab:table2}, $\theta$ and $\delta_1$ being the boundary layer momentum and displacement thickness, respectively.

\begin{table}[hbt!]
\caption{Approximate boundary layer thickness and boundary layer based Reynolds number in the central plane $z=0$.} \label{tab:table2}
\begin{tabular}{lcc}
\hline
\\[-2.7ex]
\hline
& $\Rey_\theta$ & $\delta_1$ \\
\hline
\\[-2.7ex]
\hline
\\[-2.5ex]
Before first ramp ($x=0$) & 1500 & 0.00077~\unit{m} \\
At bottom wall separation point (see Figure~\ref{fig:meanflow}) & 6875 & 0.00169~\unit{m} \\
\hline
\end{tabular}
\end{table}

\subsubsection{Reference numerical simulation setup}

In order to conduct the DNS and the legacy RANS simulations, a structured mesh is created. The characteristic sizes of the meshes are presented in Table~\ref{tab:tablemesh}, the wall units sizing are computed from an initial SA QCR RANS computation, and the maxima are given compared to the local wall unit of the closest wall, ensuring good resolution in the entire domain. Stream-wise cell sizes are adapted using linear laws, whereas wall-normal ones are using a collection of geometrical laws. The final DNS mesh contains more than $3.5$ billion points without ghost cells. The RANS mesh uses the same construction, but with much less resolution.

\begin{table}[hbt!]
\caption{Mesh characteristics.}\label{tab:tablemesh}
\begin{tabular}{lccccccc}
\hline
\\[-2.7ex]
\hline
& $n_x$ & $n_y$ & $n_z$ & $n_{\tot}$ &$\Delta x^+_{\max}$ & $\Delta y,z^+_{\max}$ & $\Delta y,z^+_{\smash{\wall}}$\\
\hline
\\[-2.7ex]
\hline
\\[-2.5ex]
DNS & $9807$ & $605$ & $605$ & $3.59\times10^9$ &$<12$ & $<10$ & $<1.4$\\
RANS & $442$ & $195$ & $195$ & $16.8\times10^6$ &$<300$ & $<100$ & $<1.4$\\
\hline
\end{tabular}
\end{table}

The DNS (see Figure~\ref{fig:introDNS}) is performed using the open-source HPC oriented FAST compressible Navier--Stokes solver from ONERA~\cite{fast}. The solver uses a finite volume approach on multi-block structured mesh using hybrid MPI/OpenMP parallel programming, cache blocking, and vectorization to improve performance. Inviscid fluxes are solved using an extended version of the scheme described in~\cite{AUSMONERA}, details about the spatial scheme are given in Appendix~\ref{app:TODO}. The viscous fluxes are computed using a second-order accurate centered scheme. Time integration is performed via an explicit low-storage third-order 3-step Runge--Kutta scheme; the time step is set to $\Delta t = 3\times10^{-9}~\unit{s}$ to ensure a CFL number lower than $0.5$ in the whole domain. The molecular viscosity is computed using Sutherland's law. The computation was run using a hybrid MPI/OpenMP parallelism on 64 computing nodes (3072 cores), using 1024 MPI processes and 3 OpenMP threads per process. The mean flow is averaged over $2\times10^6$ iterations, which corresponds to a physical time of 6~\unit{ms}.

\begin{figure}[hbt!]
\includegraphics[trim= 0cm 0cm 0cm 5cm,width=\textwidth,clip=True]{figintro.png}
\caption{Illustration of the simulation thanks to instantaneous numerical pseudo-schlieren in a constant $z$-plane ($z=0.01$) of the configuration and instantaneous stagnation pressure contour in various $x$-plane.}\label{fig:introDNS}
\end{figure}

The legacy RANS simulations are performed using the elsA software (Safran--ONERA property~\cite{elsA}) on the structured mesh presented before (see Table~\ref{tab:tablemesh}). The inviscid fluxes are computed using the AUSM+up of Liou~\cite{liou2006sequel} and an unlimited third-order MUSCL reconstruction. The viscous fluxes are solved using a centered scheme on five points, and the time integration to find the steady solution is done implicitly with a local time-stepping approach. In order to obtain the correct solution, the simulation is initialized with a subsonic region in a subpart of the diverging section and free-stream values elsewhere, the same initialisation strategy is used for the mesh adapted simulation.

\begin{figure}[h!]
\includegraphics[width=0.9\textwidth]{comp_ro.png}
\caption{Mach number distribution in the central $z$-plane. The DNS results correspond to time-averaged flow.}\label{fig:meanflow}
\end{figure}

\begin{figure}[hbt!]
\includegraphics[width=\textwidth]{initial.png}
\caption{Starting mesh for the mesh adapted RANS simulation.}\label{fig:initial}
\end{figure}

\begin{figure}[hbt!]
\includegraphics[width=\textwidth]{meshes_AIAA.png}
\caption{Meshes from the adapted (top) and elsA (bottom) SA QCR simulations in the central plane for the shock-train region.}\label{fig:STmeshes}
\end{figure}

\subsubsection{Results}

For this case, given the slow time convergence of internal flow, the tool is run for 33 adaptation steps with a varying number of time iterations, to reach a final mesh of 4~million nodes (to be compared to the 17~million cells of the legacy setup). The time-averaged Mach number fields (in the central plane) for the DNS and both RANS simulations are presented in Figure~\ref{fig:meanflow}. The flow displays all the main expected features: first, a boundary layer develops, then encounters a first shock wave as it reaches the compression ramp, followed by an expansion at the entrance of the diffuser. The flow ends up separating via a SBLI and is slowed down by a shock train to finally reach a purely subsonic state. Both RANS simulations display results that are surprisingly close to the DNS reference, given the complexity of the flow (one can note here that the use of the QCR strongly improves the quality of the RANS solution, as described in Appendix~\ref{app:turbu}). From an engineering standpoint, both methods are acceptable for vehicle design. This validates the good behavior of the toolset for turbulent internal aerodynamics. Upon closer examination of the solutions, it is evident that the elsA results show a shock train located closer to the top wall than in the mesh-adapted simulation. Consequently, the latter is in better agreement with the DNS reference. This may be due to the traditional meshing procedure used for the standard simulation, where the mesh is highly refined near the walls to properly capture the boundary layer, leading to a mesh-induced displacement of the shock train. As the refinement of the mesh-adapted simulation is fully feature-driven, no such best practice is used (see the initial mesh presented in Figure~\ref{fig:initial}), which leads to a better resolution of the mesh further away from the wall. The difference between both final meshes in the central plane for the shock train region is displayed in Figure~\ref{fig:STmeshes}. The adapted mesh clearly displays a refined region further away from the wall in addition to the near-wall refinement for the boundary layer, which may allow a better capture of the physics of the flow.

\section{Demonstration on a full vehicle}

\begin{figure}[h!]
\reflectbox{\includegraphics[width=0.7\linewidth]{japhar1_open.png}}
\caption{Skin-friction, Mach number slice in the central plane and temperature slice at the outlet illustrating the external aerodynamics of the full vehicle.}\label{fig:mock1}
\end{figure}

\begin{figure}[h!]
\includegraphics[width=0.7\linewidth]{japhar2_open.png}
\caption{Skin pressure distribution and Mach number slice in the central plane illustrating the internal aerodynamics of the full vehicle.}\label{fig:mock2}
\end{figure}

Now that the ability of the toolset to capture both complex internal and external high speed aerodynamic flowfield has been documented, a demonstration on a full mock-up vehicle can be conducted. The case consists of a high Mach number cruise of a scramjet powered vehicle. The combustion chamber is not computed but replaced by a 0D model extending from the air intake throat to the nozzle entrance. This simulation only serves as a demonstration of how we can leverage the ability of the toolset to handle complex geometry at no additional cost. Setup time for such a simulation requires only few minutes of human work. The toolset is then run on 600 Intel Cascade Lake cores for a single job of 15 hours, which allowed conducting 12 adaptation/simulation steps to reach convergence of the aero-propulsive balance, the final mesh has 35~million nodes. Samples of results are presented in Figures~\ref{fig:mock1} and~\ref{fig:mock2}. Although quantitative data cannot be presented here, the objective is to demonstrate the tool's ability to manage complex geometries (e.g., control surfaces and bleed systems). A key advantage is that geometric variability is fully transparent to the user, for example, the same simulation can be performed with different or additional control surfaces simply by modifying the input CAD file.

\section{Conclusion}

This study presents the development and validation of a fully automated, CAD-to-post toolset designed for Navier--Stokes and Reynolds-Averaged Navier--Stokes (RANS) simulations of complex hypersonic vehicles. The toolset integrates a 2D, 2D-axi and 3D vertex-centered finite volume Navier--Stokes solver, SoNICS, with an anisotropic mesh adaptation framework based on the REFINE toolbox. By leveraging automated CAD handling, surface and volume meshing, and metric-based remeshing, the toolset eliminates the need for manual intervention in the meshing process, significantly reducing setup time and human error. The workflow is entirely parallelized and managed through a Python/CGNS interface, ensuring scalability and efficiency on modern high-performance computing architectures. The toolset was rigorously validated using a series of open benchmark cases representative of the complex flow physics encountered in hypersonic applications. First, the axisymmetric cone-cylinder-flare geometry, a canonical reentry vehicle configuration, demonstrated the toolset’s ability to accurately capture shock-boundary layer interactions and flow separation, with results closely matching high-order reference solutions. Second, the BOLT forebody case, characterized by its three-dimensional geometry and intricate vortical structures, highlighted the toolset’s robustness in predicting wall heat fluxes and complex laminar boundary layers at high Reynolds numbers. The comparison with experimental and numerical data from the literature confirmed the toolset’s reliability in resolving challenging external aerodynamic features. For internal aerodynamics, a novel validation case was designed to replicate the flow conditions within a ramjet-powered vehicle’s air intake. This case included a compression ramp, diffuser, and shock-train formation, all of which are critical for engine performance evaluation. The toolset’s results were benchmarked against a Direct Numerical Simulation (DNS) and a legacy RANS simulation using the elsA solver. The comparison revealed that the mesh-adapted RANS simulations, particularly when using the Quadratic Constitutive Relation (QCR) turbulence model, provided results remarkably close to the DNS reference. This underscores the toolset’s capability to capture key flow features such as shock-boundary layer interactions, corner vortices, and shock-train evolution, even in highly turbulent and separated flow regimes. Finally, a demonstration on a full mock-up of a scramjet-powered cruise vehicle illustrated the toolset’s scalability and adaptability to complex geometries. The simulation, which included both external and internal aerodynamics, was set up with minimal human effort and efficiently converged on high-performance computing resources. This showcases the toolset’s potential for industrial applications, where rapid design iterations and high-fidelity simulations are essential. The results presented in this study validate the toolset’s ability to handle the full spectrum of hypersonic flow complexities, from external aerodynamics to internal propulsion systems. The automated workflow, combined with anisotropic mesh adaptation, ensures both accuracy and computational efficiency, making it a valuable asset for future hypersonic vehicle design and optimization. Future work will focus on extending the toolset’s capabilities to include moving bodies (through CAD manipulation), hybrid meshes, multi-species flows, advanced turbulence modeling, and GPU-accelerated remeshing to further enhance its performance and applicability. Additionally, the integration of adjoint-based optimization and data assimilation techniques will enable more sophisticated design and analysis workflows, paving the way for next-generation hypersonic vehicle development.

\section*{Acknowledgments}

The manuscript was written through contributions of all authors.

Philippe Duveau and Pierre Grenson are acknowledged for the discussion about internal aerodynamics and help to define the reference geometry, as well as Sebastien Deck, Nicolas Renard, Jaime Vacquero, Julien Husson and Lucas Manueco for the discussions on turbulence modelling and high-fidelity simulations. The authors would like to thank Bastien Andrieu for his constant effort in Paradigm improvement and for all the discussions on geometry, mesh and partitioning. Last, the authors would like to thank Ivan Mary for the help and support regarding the use of the FastS solver and numerical scheme implementation.

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

\appendix

\section{Geometry description for the internal aerodynamics case}\label{app:geom}

This appendix describes the geometry that was designed to build the internal aerodynamics reference case.

\begin{figure}[hbt!]
\includegraphics[width=1\textwidth,clip=True]{geometry.png}
\caption{3D representation of the geometry under study.}\label{fig:geom}
\end{figure}

The cross plane is top/bottom symmetrical (see Figure~\ref{fig:geom}) consisting of first a 10~\unit{cm} long flat plate to allow for the formation of a proper turbulent boundary layer (the method used to ensure a turbulent inflow is presented in Appendix~\ref{app:turbu}). Then, the internal compression is caused by a 7~\unit{cm} long, \ang{5} compression ramp. This ramp is followed by a long \ang{-2} diverging section that spans over 25~\unit{cm}. This section is representative of the diffuser in a ramjet air intake, and its angle was selected to prevent boundary layer thickening from further restricting the inviscid flow area. Finally, the flow enters a fast-diverging section with a \ang{-5}~angle, extending over 5~\unit{cm}. To ensure a well-posed boundary condition, back pressure is generated by incorporating a choked converging-diverging nozzle downstream of the diverging section. This configuration enables the use of a supersonic outflow boundary condition. The nozzle shape is defined using a Hermitian spline. An analytical (or discretized) description is available upon request. The throat closure has been chosen so that the isolator behaviour is close to critical for the flow conditions under study (terminal compression near the start of the diffuser). However, being too close to the critical point could have led to the unstart of the air intake either during DNS transient or when using some turbulence models. To alleviate this problem, the final condition under study is what could be considered an off-design condition, where the air intake does not provide the maximal efficiency. The total height of the inlet plane in 4~\unit{cm}. The 3D geometry is then created by extruding that plane over 4~\unit{cm}, creating a square shape at the inlet of the domain.

\section{Numerical method for the DNS}\label{app:TODO}

The approximate Riemann solver is the simplified AUSM(P):
\begin{equation}
F_c^{i+\frac{1}{2}} = U_{{1}}\frac{\mathbf{q}_R+\mathbf{q}_L}{2} - \Phi U_{\dis} \frac{\mathbf{q}_R-\mathbf{q}_L}{2} + P,
\end{equation}
$U_{{1}}$ and $U_{\dis}$ being the interface fluid velocity and the dissipation velocity similar to the original formulation of the scheme~\cite{AUSMONERA}, $\mathbf{q}$ is the conservative state variable vector and the subscript $\cdot_L$ (resp.\ $\cdot_R$) corresponds to left (resp.\ right) state computed thanks to a MUSCL reconstruction, $\Phi$ is the dissipation sensor which is discussed hereafter.

In order to limit the dissipation, Mary~\cite{AUSMONERA} has proposed to identify regions where cell-to-cell oscillations vanish in order to fine tune a 2nd order low dissipation scheme. The proposed oscillations sensor reads:
\begin{equation}
\Delta_i^k =
\begin{cases}
	-1	& \text{if $(\phi^k_{i+1}-\phi^k_{i})(\phi^k_{i}-\phi^k_{i-1})<0$},
\\	1	& \text{otherwise},
\end{cases}
\end{equation}
with $\phi^k_i$ the $k$-th component of the state vector of primitive variables at point $i$. Then the dissipation sensor is defined as:
\begin{equation}
\Phi^{i}= \max_k
\begin{cases}
	1	& \text{if $\sum_{l=i}^{i+1} \Delta_{l}^k \leq 0$},
\\	0	& \text{otherwise}.
\end{cases}
\end{equation}
If the solution is fully monotonous in the stencil, the scheme behaves as a fully centered scheme. If the solution is nonmonotonous, the dissipation is fully activated.

In this work, we proposed to extend the methodology to fifth-order MUSCL reconstruction. Due to the extent of the stencil the upwinding function of the scheme must be reworked. The smoothness of the solution is still evaluated by the number of slope changes in the stencil. If there are one, two, or more, the effective dissipation is set to an increased amount of the standard upwind scheme:
\begin{equation}
\Phi^{i}= \max_k
\begin{cases}
	1 			& \text{if $\sum_{l=i-1}^{i+2} \Delta_{l}^k < 0$},
\\	\alpha_2		& \text{if $\sum_{l=i-1}^{i+2} \Delta_{l}^k = 0$},
\\	\alpha_1		& \text{if $\sum_{l=i-1}^{i+2} \Delta_{l}^k = 2$},
\\	\alpha_0		& \text{if $\sum_{l=i-1}^{i+2} \Delta_{l}^k = 4$},
\end{cases}
\end{equation}
with $0<\alpha_i \leq 1$ and $\alpha_0 < \alpha_1 < \alpha_2$. This approach strongly limits the dissipation of the scheme in regions where the solution is smooth, while preserving the stability of the scheme in non-smooth regions.

For supersonic and hypersonic flows, in order to ensure the robustness of the simulation while preserving the high resolvability properties of the scheme in the smooth flow regions, the shock capturing based on the work of~\cite{sciacovelli} is used. The discontinuity sensor reads:
\begin{equation}
\kappa_i= \frac{10}{2}\bracks[\Bigg]{1-\tanh \parens[\Bigg]{2.5+10 \frac{\delta x_i}{c} \nabla \cdot \mathbf{u}}} \times \frac{(\nabla \cdot \mathbf{u})^2}{(\nabla \cdot \mathbf{u})^2+\abs{\nabla \times \mathbf{u}}^2+\epsilon} \times \abs[\Bigg]{\frac{p_{i+1}-2 p_i+p_{i-1}}{p_{i+1}+2 p_i+p_{i-1}}},
\end{equation}
with $\delta x_i$ the size of the cell in the direction of interest, $c$ the speed of sound, $p_i$ the pressure at the point $i$ and $\mathbf{u}$ the velocity vector. The sensor is used to switch from fifth to first order MUSCL state reconstruction, as following:
\begin{equation}
\mathbf{q}_{R,L}=(1-\kappa)\mathbf{q}_{R,L}^{o5}+\kappa\mathbf{q}_{R,L}^{o1}.
\end{equation}
In strong shocks, the reconstruction tends towards first order, while elsewhere it stays at fifth order, ensuring a proper stability of the numerical method in discontinuities. As the sensor is directional, a different coefficient is applied to the reconstruction in each direction. An instantaneous map of $\kappa$ for the $i$-direction in a constant $z$-plane ($z=0.01$) of the DNS is presented in Figure~\ref{fig:kappa}, showing that the first order reconstruction use is confined to shock-wave and does not reduce the order of the scheme in turbulence nor boundary layers.

For instance, the numerical scheme depicted above with tuning parameters $\alpha_0 = 0.05$, $\alpha_1 = 0.4$, $\alpha_2 = 0.6$ has been successfully used in~\cite{POULAIN2024106386}.

\begin{figure}[hbt!]
\includegraphics[width=1\textwidth]{SC.png}
\caption{Instantaneous map of $\kappa$ in the $i$-direction in a constant $z$-plane ($z=0.01$) of the configuration.}\label{fig:kappa}
\end{figure}

\section{Turbulent inflow generation for the DNS}\label{app:turbu}

One important aspect for such DNS is the need to inject turbulence at the inlet of the domain. In the present case it was decided not to use a standard turbulence generation method (such as Synthetic Eddy Method of Digital Inflow Filtering, which may both introduce nonphysical acoustic fluctuations which might be detrimental in a confined environment), instead a method inspired from recent advances in natural transition simulation~\cite{hader2018towards,lugrin907transition} is used. To create a quick, and low noise transition, streaks are first generated by alternating viscous and non-viscous wall boundary conditions patches at the inlet of the domain. The wavenumber of the created streaks is chosen to be as small as possible, the lower limit being linked to the resolution of the scheme. The patches are 4~cells wide and are repeated every 12~cells. The length of the patches is set to 40~cells. The created streaks represent a very efficient support for transitional instabilities while being small enough to not contaminate the turbulent flow with large coherent structures. Then, to trigger the transition, a random perturbation on the density ($\rho'=\rho\parens[\big]{1+0.002*A(i,j,k)}$, with $A$ a Gaussian filtered white noise, and $\rho$ the density field, see~\cite{lugrin907transition}) is injected into the boundary layer, the support of the noise injection is 130~cells in the streamwise direction and 15~cells in the wall normal direction, starting from the inlet and the wall respectively. As such, it is confined to the viscous region near the wall. This noise seeds all the instabilities of the streaks, which then quickly break down to turbulence. Thanks to that method, turbulent state is reached after around $\sim37 \delta_1$, with $\delta_1$ taken at the beginning of the convergent section, of flat plate with minimal noise emission. The flow then has $\sim90 \delta_1$ to regularize before the start of the actual geometry under study. This process is qualitatively illustrated in Figure~\ref{fig:TurbulGen}, showing an instantaneous snapshot of the creation of the streaks and their breakdown. Figure~\ref{fig:TurbuIsOk} presents the boundary layer profile and the streamwise velocity fluctuations statistics in the central $z=0$ plane just before the first ramp ($x=0$), both display a behavior typical of low-Reynolds number wall turbulence, with no spurious fluctuations coming from the turbulence generation method. The DNS profile is very close to the SA-QCR profile, except in the buffer region, where the model is known to behave poorly.

\begin{figure}[hbt!]
\includegraphics[width=\textwidth]{decl.png}
\caption{Streamwise velocity in the first cell in the wall normal direction, from an instantaneous snapshot of the DNS, illustrating the turbulence generation method, with first the inviscid patches to create the streaks, and then the quick breakdown to turbulence.}\label{fig:TurbulGen}
\end{figure}

\begin{figure}[hbt!]
\includegraphics[width=.8\textwidth]{turbu_isok_central.png}
\caption{Time-averaged streamwise velocity profile (top) and streamwise velocity fluctuation root-mean-square value (bottom) showing the good behavior of the turbulence generation. The profile is extracted from the time-averaged flow just before the compression ramp ($x=0$) in the central $z=0$ plane. DNS data are sampled one point over four.}\label{fig:TurbuIsOk}
\end{figure}

\section{Impact of the turbulence model on the results}

While outside the scope of the study, different turbulence models have been benchmarked before choosing the SA-QCR one. The benchmarked RANS turbulence models are the following: the Spalart--Allmaras~\cite{SPALARTAL} model with the compressibility correction of~\cite{DeckSA2002} (labeled SA) and the SST model of Menter~\cite{menter1992improved} (labeled as $k$-$\omega$ SST). For both models, the impact of the Quadratic Constitutive Relation (labeled QCR) \cite{spalart2000strategies} is also tested. This non-linear addition is known to drastically improve the accuracy of linear eddy viscosity turbulence models for corner flows~\cite{dandois2014improvement}.

\begin{figure}[hbt!]
\includegraphics[width=.9\textwidth]{test_ro.png}
\caption{Mach number distribution in the central $z$-plane. The DNS results correspond to time-averaged flow.}\label{fig:meanflow_central}
\end{figure}

Figure~\ref{fig:meanflow_central} presents the Mach number in the central $z$-plane for each cases (the DNS correspond to time-averaged data). On the impact of the turbulence modeling strategy, one can notice that the non-QCR simulations predict a completely different separated region topology compared to the QCR ones and the DNS. The first separation point is located much more upstream, and it is caused by a weak, poorly defined shock-wave, which may be a sign that the separation comes from strong 3D effects at the corners of the geometry. The $k$-$\omega$ simulation quickly displays a larger separated region than the SA one. The use of QCR strongly improves the behavior of the turbulence models and leads to results in the central plane that are much closer to the DNS reference, with a separation point just upstream of $x=0.15$ and a separation shock that is a well-defined oblique shock. The main differences between all the QCR simulations and the DNS is that the shock-train is a bit more confined to the top part of the diffuser and the low-speed region in the lower downstream part of the diffuser is larger. However, the results of the QCR RANS simulations are surprisingly close to the reference, demonstrating once again the accuracy improvement provided by the QCR on RANS simulations in the presence of corner flow.

\printbibliography

\end{document}