ABSTRACT. Low, ideally zero, numerical dissipation is essential for high-fidelity, scale-resolving flow simulations, as it disrupts the accurate representation of inviscid kinetic energy and entropy conservation, both fundamental to compressible turbulence. In practice, however, numerical schemes for compressible flows rely heavily on dissipation to ensure stability, creating an inherent tension between accuracy and robustness. To address this challenge, we develop a physics-compatible numerical framework that preserves both kinetic energy and entropy at the discrete level. In this talk, we present the construction of robust, non-dissipative kinetic energy and entropy preserving (KEEP) schemes, along with their application to high-Reynolds-number aircraft flows and their extension to compressible multiphase flows.

ABSTRACT. Machine learning surrogate models can produce predictions orders of magnitude faster than classical Computational Fluid Dynamics (CFD) solvers, which makes them attractive for industrial applications where many design iterations are needed. Graph Neural Networks (GNNs) are particularly well suited for CFD since they naturally handle unstructured meshes. However, their local message passing operation leads to under-reaching, limiting their ability to model the long-range interactions common in fluid dynamics. To compensate, GNNs typically require good initial conditions, incurring extra computational cost to generate them. Transformer architectures offer an alternative: their attention mechanisms enable global information exchange, but at the cost of quadratic complexity, and they lack the geometric inductive bias that GNNs possess.

We investigate hybrid architectures incorporating attention layers into GNNs to leverage the strengths of both approaches: the geometric structure awareness of GNNs and the global modeling capability of attention for long-range dependencies. We explore different configurations and demonstrate improved accuracy over GNNs alone through numerical experiments on practical datasets.

ABSTRACT. High-order methods are particularly demanding in simulations that require high accuracy. The variational reconstruction can achieve arbitrarily high-order accuracy on a compact stencil, thereby overcoming the bottleneck associated with large reconstruction stencils in the development of efficient high-order finite volume methods. However, an efficient implementation of variational reconstruction typically requires precomputing and storing large matrices, which leads to substantial memory costs. Moreover, the reconstruction procedure involve free parameters whose values significantly influence the accuracy of the resulting variational finite volume method. To address these challenges, we propose a low-storage variational reconstruction and develop a machine-learning-based optimization framework to systematically tune the reconstruction parameters.

ABSTRACT. Accelerating computational simulations remains an active area of research, pursued from multiple directions across the CFD and numerical analysis communities. In all cases, the ability to generate results more rapidly directly strengthens the role of simulations in time-critical engineering workflows. One promising approach focuses on reducing the cost of core components of the computation, particularly the repeated solutions of nonlinear algebraic systems whose linearizations yield large sparse linear systems. Computational efficiency can be improved by initializing both the nonlinear and linear solver iterations with estimates that are closer to the eventual solution, thereby reducing the number of iterations required for convergence and shortening the overall time to solution. In the present work, this acceleration is achieved by leveraging data already available in memory during a simulation to train a neural network capable of predicting updates for the nonlinear and linear solvers. These predictions are then used to initialize the solvers. A critical requirement, however, is that the machine learning components themselves must execute in less time than the savings they provide. This balance can be maintained by training a lightweight neural network to achieve sufficient, but not necessarily high, accuracy. Such models are smaller, faster to train, and still effective at improving solver convergence. The approach demonstrated here employs in-situ graph neural network training and inference to initialize both nonlinear and linear solvers. Results from one-dimensional and two-dimensional test cases show improved convergence behavior and, importantly, an overall reduction in total simulation time. The effort will expanded to investigate production-oriented solvers for the final paper.

ABSTRACT. We improve the performance of Physics-Informed Neural Networks (PINNs) for modeling discontinuous solutions in hydrodynamics and relativistic hydrodynamics. Standard PINNs often struggle with convergence and accuracy near shocks. Building on Locally-Linearized PINNs (LLPINNs), which modify the Jacobian in regions of strong compression, we remove the need for a priori shock speeds by dynamically computing them from neighboring states and enforcing entropy conditions. We also introduce Locally-Roe PINNs (LRPINNs), which employ a Roe-averaged Jacobian to enhance conservation and shock resolution. The approach is extended to two-dimensional Riemann problems using divergence-based shock detection and dimensional splitting. Compared with a high-order WENO solver, our method captures sharper shock transitions, though small-scale vortical structures remain less resolved, motivating future improvements.

ABSTRACT. Introduction

Over the last decades, due to the growing availability of computational resources, high-fidelity simulations have gained increasing interest in the prediction of compressible turbulent flows over complex geometries, such as airfoils, compressor intakes, and turbine cascades. The simultaneous presence of pressure gradients and surface curvature has a strong influence on the resulting Shock Wave-Boundary Layer Interaction (SWBLI) which can be investigated on the benchmark case of the bump. Despite its geometric simplicity, the flow over a bump has been widely studied both numerically and experimentally since it reproduces many of the key physical mechanisms observed in the SWBLI in transonic flows over supercritical airfoils. Accurately capturing these complex flow features represents one of the major challenges in the numerical simulation.

Reynolds Averaged Navier-Stokes-based methods have shown difficulties in predicting the turbulent statistics and the onset of separation, while Direct Numerical Simulation (DNS) cannot be employed at high Reynolds conditions due to the prohibitively large number of grid points required. Although Wall-Resolved Large Eddy Simulations (WRLES) are more affordable than DNS, the resolution of the inner part of the boundary layer is still computationally demanding. For this reason, Wall-Modeled Large Eddy Simulation has emerged as a suitable compromise between accuracy and computational cost, solving the outer layer as a classical LES and modeling the near-wall region.

While the large majority of CFD solvers use body-fitted grids to perform accurate simulations on complex geometries, Immersed Boundary Method (IBM) based solvers, which rely on Cartesian meshes, have shown significantly improved efficiency on modern High Performance Computing (HPC) systems. However, IBM has usually been applied in low Reynolds conditions as a result of the difficulty in clustering points in the near-wall region.

In this work, IBM-based WMLES are applied to investigate the transonic shock boundary layer interaction on a circular bump. Different simulations are performed to evaluate the capability of the current methodology in capturing the main flow physics of the phenomenon. By comparing results with experimental data, different grid resolutions are employed to effectively assess the accuracy of the approach. Finally, an analysis of the SWBLI, pressure distributions, separation, and turbulent statistics is carried out in order to identify the key physical mechanisms governing the flow at a high Reynolds number.
 

Methodology

The solver used in the present study is URANOS (Unsteady Robust All-around Navier-StOkes Solver), which integrates a dynamic wall model approach with IBM, providing an efficient and scalable framework for GPU-based architecture. The solved governing equations are the compressible Navier-Stokes equations written in conservative form, closed through the ideal gas equation of state and the constitutive expression for the internal energy.

The subgrid scale viscosity model adopted is the wall-adaptive large-eddy (WALE), while spatial discretization is performed on a structured Cartesian grid using high-order finite-difference schemes. Convective fluxes are evaluated with shock-capturing formulation, while viscous terms are discretized using central schemes. Time advancement is carried out via an explicit multi-stage Runge–Kutta method where the time-step is dynamically varied to satisfy the Courant–Friedrichs–Lewy (CFL) stability condition. In order to avoid prohibitively high computational cost related to the resolution of the inner layer, an equilibrium-based wall-stress model is used. Convective and pressure terms are assumed to be zero near the wall, so the mean wall-parallel velocity and mean wall temperature distribution are determined from the resulting ODE system. This system of equations is solved on an independent grid starting from the wall, where velocity and temperature are set by boundary conditions, to an outer interface where the value of these variables is prescribed by LES. Once the ODE system is solved, the mean wall-parallel velocity and mean wall temperature distribution are used to evaluate the wall shear stress and wall heat flux.

Regarding the surface treatment in IBM, the Ghost Point Forcing Method is employed to represent the effects of the embedded geometry. First, a properly defined mapping function classifies each grid point as fluid, ghost, or solid. For each ghost point, a corresponding image point is built in the direction normal to the body boundaries, where the flow field variables are interpolated from surrounding points. Subsequently, this value is used to enforce the wall condition by updating the value at the associated ghost point. Then, the image point can be used as the initial location of the interface between LES and wall-model allowing the simulation at high Re numbers. Further details on the algorithm adopted in the implementation of the dynamic wall-model with the IBM are shown in previous works.

Results and Discussion

Simulations were performed on a circular arc bump and then compared the results with a previous experimental study at Mach number M=0.72 and Reynolds number Re=1.6∙10⁶ based on the bump length.

The thickness of the boundary layer at the inflow is used as reference length, so computational domain is a tridimensional box with size L_xxL_yxL_z= (240x65.5x10). On the y direction the domain is extended to include the IBM block which occupies the lower part of the domain and places a bump in the center of the lower surface. The height of the bump is h = 5, and the radius is equal to 163 , with a total length of 80. The length of the up- and downstream flat plate is equal to 80. Preliminary studies highlighted the importance of an appropriate resolution in the shock and separation region. For this reason, the grid is locally refined in the near wall region spanning over the entire bump height as well as near the trailing edge in the streamwise direction. In the vertical direction, a uniform fine grid is used near the wall, transitioning to a coarser distribution in the outer region. Grid points are distributed uniformly in the spanwise direction. This choice is motivated by the necessity to ensure the correct resolution across the bump, since it is not possible to refine in normal-to-wall direction when the IBM is adopted.

The proposed IBM-WMLES approach captures the main physical phenomena of the transonic flow over a circular bump. Flow accelerates from the leading edge of the bump due to the favorable pressure gradient originating a supersonic region which terminates in the rear part of the bump with a shock wave. Close to the wall, the shock presents -shaped structure with the left-leg at x ≈ 20, whereas further from the wall it becomes nearly normal to the incoming flow direction. The pressure rise associated with the surface curvature and the interaction with the shock, determine the separation and downstream reattachment of the boundary layer. Despite the simplicity of the wall-model adopted, the analysis of mean wall pressure reveals that the shock position is predicted reasonably well in comparison with experimental data. However, the major discrepancy with respect to the reference case is in the post-shock region where equilibrium assumption of the wall-model is not satisfied.

After an initial validation with a coarse grid of approximately 37 million grid points, the accuracy and the robustness of the current approach is investigated performing simulations with two finer resolutions. Grid resolution is increased in the three dimensions doubling and triplicating the total number of points. An analysis of the turbulent statistics, pressure distribution and skin friction coefficient indicates that wall model guarantees a proper resolution of the flow field with affordable computational cost.

ABSTRACT. Introduction Aerodynamic decelerators are essential components for ensuring feasible and controlled atmospheric entry trajectories of space vehicles. Parachutes, in particular, exhibit complex dynamic behaviors during deployment and inflation, which have a strong impact on the stability of the vehicle. Despite past and ongoing efforts, correlating the performance of the parachute with the operating conditions remains a significant challenge [1]. Therefore, research efforts are still necessary to improve the understanding of parachute dynamics in light of the evolving characteristics of modern spacecraft and the wide range of atmospheric conditions associated with mission destinations. The objective of this study is to validate the recently developed numerical methods for simulating three dimensional bodies in compressible flow through benchmark cases.

Methods The computational framework is based on STREAmS 2, a Cartesian finite-difference solver for direct numerical simulation of high-speed turbulent flows [2]. The spatial discretization of the convective terms in the Navier–Stokes equations relies on a combination of a robust energy-preserving scheme and a shock-capturing scheme. The flow solver is advanced in time using a low-storage third-order Runge-Kutta scheme.

The immersed body is introduced into the framework as a triangulated (and possibly open) surface. Boundary conditions on the body are imposed at each time step using a sharp-interface Immersed Boundary (IB) method, specifically designed for thin-walled open structures. The discretization scheme is properly modified in the vicinity of the body to prevent any finite-difference stencil from containing nodes from both sides of the immersed interface. The solver relies on a multi-task implementation using Open MPI, and it is accelerated on GPUs through CUDA.

Results Figure 1 shows the pressure coefficient for a supersonic flow past circular cylinder (Ma=2, Re=300). The results agree well with those obtained in [3]. Further comparisons on both 2D and 3D test cases will be presented.

References [1] Luca Placco, Giulio Soldati, Matteo Bernardini, and Francesco Picano. On flight instabilities of capsule-rigid parachute system during supersonic planetary descent. Aerospace Science and Tech-nology, 160:110026, 2025. [2] Matteo Bernardini, Davide Modesti, Francesco Salvadore, and Sergio Pirozzoli. Streams: A high-fidelity accelerated solver for direct numerical simulation of compressible turbulent flows. Computer Physics Communications, 263:107906, 2021. [3] Daniel Canuto and Kunihiko Taira. Two-dimensional compressible viscous flow around a circular cylinder. Journal of Fluid Mechanics, 785:349–371, 2015.

ABSTRACT. To implement the k-omega shear stress transport (SST) model with the immersed boundary method, the non-plausible behaviors of the transport variables in the near-wall region are challenging. To avoid numerical problems, we adopted the inverse of omega as the transport variable and derived the model equation from the original omega equation. Additionally, we introduced a wall-damping function so that the near-wall solutions of the transport variables are linear or constant with respect to the wall distance. The derived model was implemented in the Cartesian-grid-based flow solver UTCart and validated through flow simulations around a NACA 0012 airfoil. The results of the validation problem showed good agreement with the reference results using the original k-omega SST model with body-conforming grids. The presentation will include other validation problems, e.g., the flow simulations around the full aircraft configuration, and comparisons to the results using SA turbulence model.

ABSTRACT. One of the most important and time consuming parts of a Computational Fluid Dynamics workflow is the generation of the computational mesh. Cartesian cut-cell methods circumvent this problem by automatically generating meshes without user intervention. Traditionally, cut-cell methods have struggled to provide solutions to the Reynolds-Averaged Navier-Stokes equations for high Reynolds number flows. Firstly, the irregular cell sizes near no-slip walls due to cut cells introduce noise into gradient reconstruction in the boundary layer. Secondly, satisfying the off-wall meshing requirements for wall resolved flows places too many isotropic cells in the wall-transverse direction; incurring significant cost compared to boundary-conforming methods. In this work, we compare two interior forcing point methods: the Spalart-Allmaras wall function and an ordinary differential equation (ODE) wall model. We use three NASA turbulence modeling resource verification cases: the 2D turbulent bump, the NACA 0012 airfoil, and the high-lift multielement airfoil. The ODE wall model demonstrates significant improvement over the SA wall function and comparable error convergence rates to state of the art boundary conforming methods.

ABSTRACT. We develop Global Stability Analysis (GSA) of multi-step integration methods to analyze the occurrence of physical and spurious modes arising in spectrally discretized convection-diffusion dynamics. Stability conditions are obtained and further quantified concerning the presence of both physical and spurious modes, as well as physical modes only. The Adams-Bashforth and Adams-Moulton methods provide specific examples for which GSA yields explicit conditions on the Peclet and CFL numbers that need to be respected to achieve simulations free of spurious modes.

ABSTRACT. Convection-diffusion equations constitute the mathematical foundation of many governing laws in Computational Fluid Dynamics (CFD), including the Navier-Stokes and Reynolds-Averaged Navier-Stokes (RANS) equations. To numerically resolve these equations, High-Order (HO) schemes have gained significant popularity in recent years due to their computational efficiency in resolving smooth regions of the solution.

Despite these advantages, achieving robust HO discretizations remains challenging. While stabilization for hyperbolic operators is well established, robust treatment of diffusion—especially for strongly anisotropic diffusion on distorted meshes—remains an active research area. Moreover, the parabolic nature of the problem imposes severe time-step restrictions on explicit schemes due to the von Neumann stability condition, making fully implicit formulations essential.

In this work, we present a theoretical analysis of a fully discrete, time-implicit DGSEM based on the Bassi and Rebay 2 (BR2) formulation, focusing on domain invariance and entropy stability. We derive robustness conditions and introduce an algebraic graph viscosity for stabilization. A Flux-Corrected Transport (FCT) limiter is used to couple a robust low-order scheme with the high-order operator, preserving accuracy while enhancing stability. These properties are extended to generic distorted meshes, and numerical benchmark tests are provided to validate the theoretical results.

ABSTRACT. This paper addresses the need for time-continuous temporal predictors in time-accurate local time stepping (LTS) for high-order discontinuous Galerkin (DG) methods. Within a predictor--corrector DG--LTS framework, we propose a family of multistep, multistage continuous explicit Runge--Kutta predictors (MSCERK). The key idea is to incorporate information from previous steps into the continuous-extension structure of the current step, distributing the information required for a time-continuous representation across multiple steps. For a prescribed target order, MSCERK reduces the minimum number of stages that must be explicitly constructed in each step to achieve the same continuous order, thereby decreasing predictor-stage right-hand-side (spatial-operator) evaluations and related local costs in LTS settings while retaining the structural benefits of time-continuous representations for asynchronous interface-flux integration and interpolation. A linearly stable, third-order accurate two-step two-stage scheme, MSCERK(2,2), is presented as an example. Numerical results demonstrate approximately fourth-order convergence for a \(P^3\) Burgers' equation benchmark and show that MSCERK(2,2) reproduces reference wall-pressure spectra in the \(30p30n\) benchmark case while reducing the computational cost by 27.1\% with negligible changes （less than 0.7%）in integral aerodynamic coefficients.

ABSTRACT. 1. Introduction Developments toward entropy-stable high-order methods have brought high-fidelity computational fluid dynamics (CFD) simulations closer to industrial application. In particular, entropy-stable discretizations have demonstrated exceptional robustness for simulations of unsteady flows, even in the presence of complex geometries. Nonlinear stability properties enabled by entropy stability – first developed by Harten [1] – and discrete summation-by-parts properties [2] have demonstrated robustness properties superior to standard high-order discretizations of conservative forms [3, 4, 5]. While high-order nonlinearly stable discretizations offer robustness advantages, the best choice of timestepping scheme remains unclear. Explicit Runge-Kutta (RK) is a common choice; others use implicit-explicit (IMEX) [6] or a variety of implicit methods, for instance in [7]. In this work, we will compare temporal integration methods available in the Parallel High-Order Library for PDEs (PHiLiP), an open-source code developed by the Computational Aerodynamics Group at McGill [8]. We will aim to compare the properties of explicit RK, diagonally-implicit RK, and a space-time discretization on an airfoil test case. 2. Methodology We use the nonlinearly stable flux reconstruction (NSFR) scheme of Cicchino, Nadarajah et al.[9, 10, 11] as a foundation for the results of this abstract. We direct readers to the papers [9, 10, 11, 12] for a detailed description of the method, and only highlight key features in the interest of brevity. The NSFR scheme is arbitrarily high-order and benefits from larger time steps due to flux reconstruction, which is applied as a modal filter according to [13, 14]. The skew-symmetric stiffness operator of Chan [15] enables entropy stability while allowing a wide range of node choices. We use high-order, curvilinear elements following [16]. The implementation scales at $O(p^{dim+1}$ [17] because of the use of sum-favctorized tensor product elements [18]. The PHiLiP code code also supports the use of a strong form DG spatial discretization, which uses similar sum factorization routines. The two options are similar in code structure and allow for even comparison. 2.1 Explicit Runge-Kutta Temporal Integration The typical choice of time stepping method for high-order methods is an explicit $s$-stage RK method, which is defined according to a Butcher tableau with coefficients $a_{ij}$ and weights $b_i$ [19]. The relaxation RK method of Ketcheson [20] and Ranocha [21] can be used to ensure entropy stability. Explicit methods have a strictly lower-triangular coefficient matrix. While RK methods offer convenient implementation, the high aspect ratio cells found in wall-bounded flows lead to very small stable time step sizes, particularly considering that the CFL number itself is inversely proportional to the grid’s polynomial order [22]. 2.2 Diagonally-Implicit Runge-Kutta Temporal Integration A wide variety of implicit RK methods can be used to circumvent the small timesteps required by wall-bounded flows. Here, we focus on diagonally-implicit RK, which has nonzero $a_{ij}$ on and below the diagonal, allowing each stage to be solved sequentially [23]. Presently, we use a Jacobian-free Newton-Krylov method without preconditioning for implicit solves [24, 25]. Convergence tolerances are chosen such that relatively few GMRES iterations are used per nonlinear iteration, which we find to reduce cost per time step. 2.3 Space-Time Nonlinearly Stable Flux Reconstruction The third option we investigate for discretizing the unsteady term is a space-time method, where time is treated as an additional dimension, and the resulting system is solved implicitly. Space-time elements use the same polynomial bases for spatial dimensions and the added temporal dimension. The major benefit of space-time discretizations is the ease of $hp$ adaptation in the temporal dimension: temporal order and number of timesteps can be locally adapted using DG mechanisms. For instance, the work of Zahr [26] has used $r$ adaptation in space-time to track moving shocks. However, implementing space-time into an existing code can be prohibitively difficult. Thus, the present work will be a significant development toward the understanding of how space-time discretizations compare to more conventional choices. The space-time method used in this work is inspired by that of Friedrich et al. [27]. Similar skew-symmetric forms are used in both space and time, enabling fully-discrete entropy stability and preservation proofs. Likewise, we also formulate a strong DG space-time discretization. The method is formulated to use purely upwind numerical flux at interfaces between time-slabs, allowing the implicit solve to treat slabs in a decoupled fashion. The order of accuracy and entropy stability properties have been verified in the PHiLiP implementation, but airfoil results are not included in this abstract. 2.4 NACA0012 Test Case We define a test case which allows for the comparison of the temporal schemes introduced previously. We seek a test case which mimics a wall-resolved turbulent flow, where unsteadiness must be captured in a robust manner, and cells must be very small in the near-wall region. We thus define our test case as flow around a NACA0012 airfoil at $Re=1E5$ and $8$-degree angle of attack, resulting in an inherently unsteady flow. The flow characteristics are $M_{\infty} = 0.1, \gamma= 1.4, Pr= 0.72, T_{\infty} = 273.15$, corresponding to nearly-incompressible flow. The case is run in two dimensions because of the grid limitations for the spacetime case. The grid is very fine in the near-wall region, simulating the setup one would use in wall-resolved large eddy simulation. We aim to use this test case to compare the computational cost of explicit and implicit time stepping schemes while retaining consistent flow characteristics. Furthermore, the inclusion of strong DG and NSFR allows for comparisons between entropy-stable and conventional high-order methods in the context of implicit methods. 3. Conclusions This work will compare timings using explicit RK, implicit RK and space-time temporal discretizations using an unsteady airfoil test case. The work will compare computational cost on a consistent, well-scaling code. The test case aims to replicate the computational challenges of wall-resolved LES and highlights a case where implicit methods may be needed. We will assess the properties of the three options above and give a cost-based recommendation for further development. 4. Acknowledgements We acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant Program [RGPIN-2019-04791] and Canadian Graduate Scholarships - Doc- toral program [CGS-D-579552], the support of McGill University, and the Zonta International Amelia Earhart Fellowship program. C.P. thanks Dominic Roy and Harsh Kumar for their aid with the mesh and test case setup. References [1] Ami Harten, Bjorn Engquist, Stanley Osher, and Sukumar R Chakravarthy. Uniformly high order accurate essentially non-oscillatory schemes, III. Journal of Computational Physics, 71(2):231–303, 1987. [2] Travis C Fisher, Mark H Carpenter, Jan Nordström, Nail K Yamaleev, and Charles Swanson. Discretely conservative finite-difference formulations for nonlinear conservation laws in split form: Theory and boundary conditions. Journal of Computational Physics, 234:353–375, 2013. [3] Gregor J Gassner, Andrew R Winters, and David A Kopriva. Split form nodal discontinuous galerkin schemes with summation-by-parts property for the compressible euler equations. Journal of Computational Physics, 327:39–66, 2016. [4] Diego Rojas, Radouan Boukharfane, Lisandro Dalcin, David C Del Rey Fernández, Hendrik Ra- nocha, David E Keyes, and Matteo Parsani. On the robustness and performance of entropy stable collocated discontinuous galerkin methods. Journal of Computational Physics, 426:109891, 2021. [5] Jesse Chan, Hendrik Ranocha, Andrés M Rueda-Ramírez, Gregor Gassner, and Tim Warburton. On the entropy projection and the robustness of high order entropy stable discontinuous Galerkin schemes for under-resolved flows. Frontiers in Physics, page 356, 2022. [6] Alex Kanevsky, Mark H Carpenter, David Gottlieb, and Jan S Hesthaven. Application of implicit– explicit high order Runge–Kutta methods to discontinuous-Galerkin schemes. Journal of Compu- tational Physics, 225(2):1753–1781, 2007. [7] Ralf Hartmann, Francesco Bassi, Igor Bosnyakov, Lorenzo Botti, Alessandro Colombo, Andrea Crivellini, Matteo Franciolini, Tobias Leicht, Emeric Martin, Francescocarlo Massa, et al. Implicit methods. In TILDA: Towards Industrial LES/DNS in Aeronautics: Paving the Way for Future Accurate CFD-Results of the H2020 Research Project TILDA, Funded by the European Union, 2015-2018, pages 11–59. Springer, 2021. [8] Doug Shi-Dong and Siva Nadarajah. Full-space approach to aerodynamic shape optimization. Computers & Fluids, 218:104843, 2021. [9] Alexander Cicchino, Siva Nadarajah, and David C Del Rey Fernández. Nonlinearly stable flux Thirteenth International Conference on Computational Fluid Dynamics (ICCFD13) Milan, Italy, July 06-10, 2026 4 reconstruction high-order methods in split form. Journal of Computational Physics, 458:111094, 2022. [10] Alexander Cicchino, David C Del Rey Fernández, Siva Nadarajah, Jesse Chan, and Mark H Car- penter. Provably stable flux reconstruction high-order methods on curvilinear elements. Journal of Computational Physics, 463:111259, 2022. [11] Alexander Cicchino and Siva Nadarajah. Discretely nonlinearly stable weight-adjusted flux recon- struction high-order method for compressible flows on curvilinear grids. Journal of Computational Physics, 521:113532, 2025. [12] Carolyn MV Pethrick and Siva Nadarajah. Fully-discrete nonlinearly-stable flux reconstruction methods for compressible flows. Journal of Computational Physics, 533:113984, 2025. [13] Y Allaneau and Antony Jameson. Connections between the filtered discontinuous Galerkin method and the flux reconstruction approach to high order discretizations. Computer Methods in Applied Mechanics and Engineering, 200(49-52):3628–3636, 2011. [14] Philip Zwanenburg and Siva Nadarajah. Equivalence between the energy stable flux reconstruction and filtered discontinuous Galerkin schemes. Journal of Computational Physics, 306:343–369, 2016. [15] Jesse Chan. Skew-symmetric entropy stable modal discontinuous galerkin formulations. Journal of Scientific Computing, 81(1):459–485, 2019. [16] David A Kopriva. Metric identities and the discontinuous spectral element method on curvilinear meshes. Journal of Scientific Computing, 26:301–327, 2006. [17] Alexander Cicchino and Siva Nadarajah. Scalable evaluation of hadamard products with tensor product basis for entropy-stable high-order methods. Journal of Computational Physics, page 113134, 2024. [18] Steven A Orszag. Spectral methods for problems in complex geometrics. In Numerical methods for partial differential equations, pages 273–305. Elsevier, 1979. [19] John Charles Butcher. Numerical methods for ordinary differential equations. John Wiley & Sons, 2016. [20] David I Ketcheson. Relaxation Runge-Kutta methods: Conservation and stability for inner-product norms. SIAM Journal on Numerical Analysis, 57(6):2850–2870, 2019. [21] Hendrik Ranocha, Mohammed Sayyari, Lisandro Dalcin, Matteo Parsani, and David I Ketcheson. Relaxation Runge-Kutta methods: Fully discrete explicit entropy-stable schemes for the compress- ible Euler and Navier-Stokes equations. SIAM Journal on Scientific Computing, 42(2):A612–A638, 2020. [22] Bernardo Cockburn and Chi-Wang Shu. Runge–Kutta discontinuous Galerkin methods for convection-dominated problems. Journal of Scientific Computing, 16:173–261, 2001. [23] Christopher A Kennedy and Mark H Carpenter. Diagonally implicit Runge–Kutta methods for ordinary differential equations. A review. Technical report, 2016. [24] Dana A Knoll and David E Keyes. Jacobian-free Newton–Krylov methods: a survey of approaches and applications. Journal of Computational Physics, 193(2):357–397, 2004. [25] Eduardo Jourdan de Araujo Jorge Filho and Zhi J Wang. A matrix-free gmres algorithm on gpu clusters for implicit large eddy simulation. In AIAA Scitech 2021 Forum, page 1837, 2021. [26] Charles J Naudet and Matthew J Zahr. A space-time high-order implicit shock tracking method for shock-dominated unsteady flows. Journal of Computational Physics, 501:112792, 2024. [27] Lucas Friedrich, Gero Schnücke, Andrew R Winters, David C Del Rey Fernández, Gregor J Gassner, and Mark H Carpenter. Entropy stable space–time discontinuous Galerkin schemes with summation- by-parts property for hyperbolic conservation laws. Journal of Scientific Computing, 80(1):175–222, 2019.

ABSTRACT. 1. Introduction The predictive simulation of reacting flows is essential for designing next-generation aerospace propulsion systems, including scramjets, rotating detonation engines, and hypersonic propulsion concepts. These applications involve strong coupling between compressible fluid dynamics, heat release, molecular transport, and finite-rate chemical kinetics, leading to thin reaction zones, sharp thermochemical gradients, shock waves and multiscale dynamics. Accurate numerical resolution is critical for predicting flame dynamics under extreme operating conditions. Yet, conventional low-order discretization methods often suffer from excessive numerical dissipation and dispersion that limiting predictive accuracy for predicting the relevant combustion regimes.

High-order discretization methods offer a promising pathway to overcome these limitations by providing superior accuracy per degree of freedom and reduced numerical dissipation, enabling improved resolution of steep gradients and fine-scale thermochemical structures. In particular, p-refinement achieves exponential convergence and lower cost compared to h-refinement, driving increasing adoption of high-order approaches. Among these, the Spectral Element Method (SEM) combines high accuracy with computational efficiency. Its element-local formulation yields favorable properties, including diagonal mass matrices and high arithmetic intensity, making it well suited for accelerated architectures based on graphics processing units (GPUs).

GPUs have emerged as a key enabling technology, providing massive parallelism and high memory bandwidth well suited for compute-intensive, element-local operations and stiff chemistry integration, which typically dominates runtime. High-order spectral element methods are particularly well aligned with GPU architectures due to their localized data access and high arithmetic intensity, enabling efficient acceleration. Efficient GPU-based integration of finite-rate chemical kinetics is also essential to fully exploit modern heterogeneous architectures. Leveraging GPUs is therefore essential to achieve the performance and scalability required for predictive high-order combustion simulations. These challenges are particularly evident across combustion regimes. Deflagrations involve thin reaction zones governed by diffusive and reactive transport, demanding high-order accuracy to resolve flame structure and propagation. In contrast, detonations involve tightly coupled shock–reaction interactions, producing strong discontinuities and complex wave dynamics that impose stringent requirements on stability, shock resolution, and computational efficiency.

While several GPU-accelerated solvers exist for discontinuous Galerkin methods and incompressible flows, GPU-native high-order continuous Galerkin solvers capable of handling compressible reacting flows with finite-rate chemistry remain largely unexplored. The present work addresses this gap by extending the GPU-accelerated spectral element solver SOD2D to simulate compressible multi-component reactive flows with finite-rate chemistry. The approach combines high-order spatial discretization with entropy-viscosity shock capturing and projection-based stabilization to ensure numerical robustness while preserving accuracy. Chemical kinetics integration is fully offloaded to GPUs through a dedicated acceleration library, enabling efficient treatment of stiff reaction mechanisms. The resulting framework enables accurate and scalable simulations of reactive flows across combustion regimes, from laminar premixed flames to detonation waves, providing a foundation for predictive simulations of high-speed propulsion systems on modern GPU-accelerated supercomputers.

2. Methodology 2.1 Numerical Methods The present work employs SEM as the spatial discretization framework. A key advantage lies in its computational efficiency: unlike the sparse banded mass matrix in traditional FEM, SEM's mass matrix is purely diagonal, significantly reducing memory requirements and computational cost. This property arises from the use of Lagrange basis functions, with interpolation and integration collocated at Gauss-Lobatto-Legendre (GLL) nodes. Since the basis functions satisfy the Kronecker delta property, for the weak form of the mass matrix, only the diagonal terms with i=j yields non-zero contributions, thus leading to a diagonal mass matrix.

The implementation employs a hybrid approach combining CUDA and OpenACC directives, balancing performance with portability and maintainability. The element-local nature of SEM naturally lends itself to parallel execution: nodal operations require only data from the corresponding element in most situation, thus minimizing global memory traffic.

The chemistry integration, which dominates cost due to stiff ODE solving, is fully offloaded to GPUs through ChemInt—a newly developed in-house library for GPU-accelerated stiff chemical kinetics. All the development and performance benchmarking above are conducted on the MareNostrum 5 at Barcelona Supercomputing Center, utilizing NVIDIA H100 GPUs for multi-GPU strong scaling.

2.2 Governing Equations and Stabilization Techniques The governing equations for compressible multi-component reactive flows comprise the Navier-Stokes equations along with the species transport equations for every species s, where t denotes the physical time, \rho, \boldsymbol{u}, p, T, \boldsymbol{\tau}, \lambda, \boldsymbol{V_k}, Y_s, \boldsymbol{V_s} and \dot{\omega}_s represent the density, velocity vector, pressure, temperature, viscous stress tensor, thermal conductivity, correction velocity, and the mass fractions, diffusion velocities, and reaction rates of species s, respectively. Additional entropy viscosity (EV) terms are added on the right-hand sides of the governing equations. These artificial diffusion terms rely on an additional viscosity \mu_\mathrm{EV}, which is constructed following the entropy residual approach described in the work of -L. Guermond, et al. Therefore, these EV terms capture shocks and only be activated nearby, thus stabilizing the flow precisely. This local low-order stabilization is completed with global high-order "Local Projection Stabilization" (LPS) terms, included following the procedure described in, when computing the equations weak form integrals within the continuous Galerkin spectral element method framework. This projection-based stabilization scales with O (h^{p+1}), becoming increasingly targeted as polynomial order increases and thus minimizing interference with resolved flow features.

Finally, the system is closed by thermodynamic relations and temperature-dependent transport properties for an ideal gas mixture. To further enhance numerical stability, a flux splitting strategy is applied to the convective terms following the work of Kennedy and Gruber. The time integration of the convective and diffusive fluxes is performed using a third-order strong stability-preserving explicit Runge-Kutta schemes, while operator splitting is applied to integrate the chemical reactions rates \dot{\omega}_s.

3. Results 3.1 One-dimensional Hydrogen Laminar Premixed Flame First, a one-dimensional laminar premixed hydrogen flame is simulated for a stoichiometric H2/air mixture, at initial conditions of 300K and 1atm, and with the reduced mechanism of P. Boivin et al. involving 9 chemical species and 12 reactions. The computational domain has been discretized with third-order spectral elements of length h_e=50~\mathrm{\mu} m, while the time step has been computed from Courant-Friedrich (CFL) and Fourier (Fo) numbers set to 0.9 and 0.3, respectively. The computed flame structure and propagation speed show excellent agreement with reference solutions obtained from Cantera (seen in Figure 1), validating the accurate implementation of transport properties and chemical kinetics.

3.2 Two-dimensional Hydrogen/Oxygen/Argon Detonation Second, to evaluate the solver's performance in capturing supersonic combustion regimes, a two-dimensional H2/O2/Ar detonation at low pressure is simulated with the same reduced mechanism from P. Boivin et al. The space-time discretization is here set-up with third-order quadratic spectral elements of size h_e=75~\mathrm{\mu} m, as well as (CFL,Fo)=(0.35,0.1) conditions. The numerical simulation is initialized following the work of R. F. Johnson et al., allowing the development of a self-sustained detonation wave. The resulting cellular detonation structure (seen in the Figure 2), characterized by triple points and transverse waves, is successfully captured, demonstrating the solver's ability to simulate the detonation initiation and quenching processes, and the ability to resolve complex shock-chemistry interactions inherent in hydrogen detonations.

4. Conclusions A GPU-accelerated reactive flow solver has been developed by extending the spectral element framework SOD2D with finite-rate chemistry capabilities. Advanced stabilization techniques ensure numerical robustness for supersonic flows with strong discontinuities while preserving high-order accuracy. Validation cases demonstrate accurate resolution of both laminar flame propagation and detonation structures, establishing a promising computational framework for predictive simulations of high-speed propulsion systems on modern GPU-accelerated supercomputers.

ABSTRACT. Flame control strategies can be broadly classified into passive and active approaches. Passive methods are simple and require no external energy input, however, they lack adaptability to changing flow conditions. Typically, they involve geometric modifications of the injection system to stabilize the flame by altering the flow field.

Active flow control methods offer greater flexibility, as they do not rely solely on optimizing the flow domain for a specific operating regime. Their primary advantage lies in their ability to dynamically adjust control parameters, such as excitation amplitude and frequency, enabling real-time responses to variations in flow conditions, for example, changes in inlet velocity or temperature.

In the present work, we employ both passive and active flow control strategies to improve the combustion characteristics of nitrogen-diluted hydrogen and ammonia–hydrogen flames, which are considered eco-friendly, carbon-free energy carriers. Injection systems incorporating cylindrical and star-shaped bluff bodies are investigated, with harmonic forcing applied to exploit their combined effects and enhance fuel–oxidizer mixing. The study is conducted using Large Eddy Simulation (LES) and an advanced numerical framework based on a two-stage procedure that involves two computational codes, each with second-order finite-volume and sixth-order compact difference schemes. The primary objective is to assess how the bluff-body shape and excitation parameters (frequency and amplitude) influence the formation of large- and small-scale vortical structures, and how the resulting flow modifications affect flame characteristics, including flame shape, fuel consumption, temperature distribution, and NOx formation.

ABSTRACT. The growing recognition of the importance of detailed reaction mechanisms in high-speed combustion flows has made their application in numerical simulations increasingly essential. To address the challenges associated with implementing these detailed mechanisms, this paper proposes a three-tier collaborative computational framework. This framework dynamically reduces the mechanism scale during computation while maintaining high-performance parallel processing. In a homogeneous reactor model, the parallel implementation of this framework achieved a 30% reduction in computational time. Furthermore, when applied to an oblique detonation flow field, the framework reduced the reaction computation time by 93%.

ABSTRACT. Initiated chemical vapour deposition (iCVD) is an increasingly important technique for the fabrication of conformal polymer thin films, offering a solvent-free process with excellent compositional control and compatibility with temperature-sensitive substrates. Despite these advantages, the iCVD reactor performance is highly sensitive to operating conditions, and process optimisation is often carried out through extensive experimental trial-and-error. In this work, a multiphysics 3-D model of a laboratory-scale iCVD reactor is developed to provide mechanistic insight into the coupled transport and reaction phenomena governing polymer film growth. The model combines laminar non-isothermal flow, multicomponent mass transport, electric potential, and detailed surface polymerisation reactions. The gas-phase initiator decomposition, precursor transport, and surface adsorption, initiation, propagation, and termination reactions are explicitly accounted for. Simulations allow to determine the spatial distributions of velocity, temperature, and chemical concentrations under low-pressure operating conditions, typical of these reactors. The radiative heat transfer from the filament is shown to dominate the thermal field, strongly influencing initiator decomposition and surface reaction rates. Predicted deposition rate profiles exhibit near-radial symmetry and high uniformity across the substrate, with minor variations of approximately 2 nm/min. A systematic comparison between simulated and experimental deposition rates over a range of pressures demonstrates good agreement, validating the modelling framework. A sensitivity analysis indicates that the surface propagation rate constant is the dominant kinetic parameter affecting deposition rates, while reactor pressure, filament voltage, cooling stage temperature, and inlet flow rates significantly influence substrate temperature. All in all, the model provides valuable insight into iCVD reactor behaviour and offers a robust tool for guiding reactor design and preliminary process optimisation, reducing reliance on the costly experimental iterations.

ABSTRACT. Finite element-based high-order methods have emerged as a powerful alternative to overcome the weakness of finite volume methods (FVM), yet they remain vulnerable to shock-induced instabilities such as subcell Gibbs–Wilbraham (G–W) oscillations. To eliminate these oscillations, shock-capturing methods have been developed by combining a troubled-cell indicator with a local stabilization, such as a limiting strategy or an artificial viscosity. Nonetheless, conventional methods are not free from a fundamental trade-off between accuracy and robustness, largely due to the coupled nature of (i) where to flag troubled-cells and (ii) how to prescribe the form and strength of the local stabilization. Recently, the shock-capturing PID controller (SPID) achieved excellent performance up to the range of moderately strong shocks via the linear control theory with the data-driven gain optimization, while its robustness and convergence significantly degrade for hypersonic flows due to the strong non-linear mapping from the PID error signal to the required artificial viscosity. To extend the capability of SPID to hypersonic flows, we develop a fully data-driven shock-capturing framework comprising (i) a parameter-free, decision tree-based troubled-cell indicator and (ii) an ANN-assisted dilatation-based artificial viscosity (DB) for strong shocks. By introducing the data-driven DB model into SPID, a hybrid artificial viscosity method (DB_SPID) is proposed, which significantly improves the performance of SPID in hypersonic flows with a target-oriented flow stabilization over a wide range of discontinuity types and strengths. Through numerical experiments for hypersonic inviscid and viscous flows, it is demonstrated that the proposed DB_SPID method significantly improves robustness and convergence characteristics in high Mach number flows.

ABSTRACT. Numerical simulations of hypersonic flows are challenging due to the combination of thermophysical modeling and numerical strategy. Specifically, in the context of atmospheric vehicle entry, the solver has to deal with strong shock waves which provoke energy excitation and molecular dissociation. This leads to a tremendous heat transfer on the vehicle surface, mainly due to conduction. However, in the presence of gas-surface interactions, mass diffusion becomes a relevant contribution to the total heat flux due to catalysis and ablation, the latter rising from the interaction between the chemical species in the boundary layer and the material composing the Thermal Protection System (TPS). Nitridation and/or oxidation processes may take place, pushing the gas far from the wall (blowing effect) and mitigating the heat flux on the surface by material decomposition and mass loss. In some case, the surface recesses with a finite velocity, provoking a shape change of the TPS.

Most of the numerical solvers in the scientific community employ body-fitted (BF) grids, which allow the use of large cell aspect ratios, thus ensuring very high wall normal resolution in the boundary layer. However, when dealing with complex/moving geometries, multiblock structured grids or unstructured grids are necessary. As an alternative, immersed-boundary (IB) techniques exploit Cartesian meshes for the fluid domain, where the geometry is embedded as a Lagrangian (solid) domain. The governing equations are then solved in the Eulerian (fluid) domain with a forcing term which accounts for the forces exchanged between the fluid and the body, such that the equations obey to the desired boundary condition on the surface (e.g., no-slip). Specifically, thanks to a ray-tracing procedure, the cells close to the surface are tag as interface cells: in these cells, an appropriate reconstruction of the flow field variables is needed to satisfy the wall boundary condition.

This work presents numerical results for hypersonic flows obtained with two different IB techniques, compared to reference results to provide an assessment about the accuracy of IB techniques under different conditions.

ABSTRACT. Aleatoric uncertainties are propagated through a computational model representative of a reference experimental hypersonic flow, which exhibits strong nonequilibrium thermochemical effects. Namely, a hypersonic stream over a double-cone axisymmetric geometry for which observations are available. The experiment is reproduced numerically. The inflow uncertainty is first characterized, and then sampled to generate a database. The set of realizations is then postprocessed to build standard Polynomial Chaos Expansions to perform the uncertainty analysis. A global sensitvity analysis is carried out through Sobol's method. Then, the true experimental operating conditions are inferred based on the available observations through a calibration procedure.

ABSTRACT. 1. Introduction The extreme temperatures generated near the surface of a vehicle moving at hypersonic speeds give rise to complex high-enthalpy processes, including vibrational excitation, chemical reactions, thermal relaxation and gas-surface interactions. In fully-developed high Mach number turbulent boundary layers, the coupling of velocity and temperature fields, combined with shock wave interactions and finite-rate chemistry effects, changes turbulence structure and transport properties [1]. For this reason, there is a direct impact on the prediction of wall quantities that are important for the design of a Thermal Protection System (TPS) and for reliable assessment of aerodynamic heating in hypersonic applications. In this context, Direct Numerical Simulation (DNS) represents the most accurate way to analyze turbulent flows, as it resolves all turbulent scales without relying on turbulence models. However, its computational cost restricts applications to canonical configurations at moderate Reynolds numbers, such as turbulent channel flows.In recent years, extensive DNS and WRLES databases for supersonic and hypersonic boundary layers have been generated, spanning a wide range of Mach numbers, wall temperature ratios and thermochemical conditions [2,3]. These databases have proved to be of fundamental importance because they provide key reference parameters for the development of wall models for LES. Recent developments have extended ODE-based wall models to high-enthalpy turbulent boundary layers, solving simplified momentum and energy equations in ordinary differential form and incorporating enthalpy–velocity relations and generalized Reynolds analogy [4]. These approaches have shown promising results for hypersonic WMLES configurations under thermochemical non-equilibrium conditions. The objective of the present study is twofold. First, to generate a high-fidelity DNS database of compressible turbulent channel flows at moderate Mach numbers, both in perfect- and reactive-gas conditions, which will serve as a consistent reference framework for hypersonic wall-model development. Second, to perform a priori and a posteriori studies of wall-model assumptions based on the generated DNS database. DNS-resolved near-wall quantities are used to directly test the equilibrium hypotheses underlying reduced wall-model formulations. The same wall-model equations are then integrated within the numerical framework and their performance is evaluated in a fully coupled simulation environment.

2. Methodology The governing equations correspond to the compressible multi-species Navier–Stokes system with finite-rate chemistry and vibrational relaxation. In conservative form, the formulation includes transport equations for each chemical species, for momentum, for total energy and for vibro-electonic energy. The species equations describe the temporal variation of the partial densities as a balance between convection, molecular diffusion and chemical production rates. The momentum equation accounts for convective transport, pressure forces and viscous stresses. The total energy equation accounts for convective and pressure work terms, viscous dissipation and heat fluxes. An additional equation governs the vibro-electronic energy mode, including its transport and source terms due to thermochemical relaxation. Spatial discretization is performed using a high-order finite-difference formulation in generalized curvilinear coordinates with one-dimensional wall-normal stretching. The convective terms are discretized in a sixth-order skew-symmetric formulation in order to preserve the discrete kinetic energy balance and reduce aliasing errors in turbulent regimes. In the presence of strong gradients, a hybrid strategy is adopted in which the skew-symmetric central discretization is locally replaced by a fifth-order WENO reconstruction to ensure robustness and shock-capturing capability. Viscous terms are discretized using sixth-order central differences. Time integration relies on a third-order accurate, explicit Runge–Kutta scheme combined with operator splitting for stiff thermochemical source terms [5]. The solver used, DAMISO, is implemented within a multi-GPU MPI-CUDA framework and has been previously verified on canonical compressible benchmarks.

3. DNS database generation The baseline DNS database includes turbulent channel flows under both low- and high-enthalpy conditions. All simulations are performed in a canonical plane channel configuration with periodic boundary conditions imposed in the streamwise and spanwise directions and no-slip isothermal walls in the wall-normal direction. A one-dimensional grid stretching is applied toward both walls to achieve a target near-wall resolution of y+ ≈ 0.8. In order to speed up convergence, simulations are first initialized on a coarser grid and subsequently interpolated onto a finer one with final DNS resolution. Time averages are collected when the bulk temperature and wall density in the domain reach a statistically steady state. Two perfect-gas configurations are considered at Mach 1.5 and Mach 3, corresponding to well-established reference results by Modesti & Pirozzoli [6]. The computational domain extends over Lx × Ly × Lz = 4πh × 2h × (4/3)πh, where h denotes the channel half-height. For the Mach 1.5 case (PG-M1.5), the friction Reynolds number is Reτ = 219.7 and the friction Mach number is Maτ = 0.079, with a heat-flux parameter Bq = 0.048. The grid resolution is 256 × 128 × 128 in the streamwise, wall-normal and spanwise directions, respectively. The near-wall resolution is y+wall ≈ 0.8, with streamwise and spanwise spacings of Δx+ ≈ 10.8 and Δz+ ≈ 7.2. For the Mach 3 case (PG-M3), the friction Reynolds number is Reτ = 448.0 and the friction Mach number is Maτ = 0.11, with Bq = 0.14. The grid resolution is 570 × 260 × 250, maintaining y+wall ≈ 0.8, with Δx+ ≈ 10.0 and Δz+ ≈ 7.6. Instantaneous fields of density-gradient magnitude for the PG-M3 case show intense small-scale structures and elongated near-wall streaks, highlighting the presence of fully developed compressible turbulence. Reynolds-averaged wall-normal profiles for velocity, temperature, viscosity and turbulence intensities display excellent agreement with literature data [6], further verifying the numerical setup. The DNS database will be enriched with two high-enthalpy cases, at Mach 3 and with wall temperature values of 1750 K and 3500 K, respectively. For these cases, the flow is initialized with constant pressure, temperature and air composition (i.e., XN2 = 0.79 and XO2 = 0.21).

4. A-priori and A-posteriori assessment of wall modeling The second objective of the present study is the assessment of wall-model assumptions for hypersonic applications. A one-dimensional equilibrium ODE-based wall model is being implemented within the numerical framework to provide a consistent bridge between the DNS database and future WMLES applications in hypersonic regimes. Under local equilibrium and constant total-shear-stress assumptions, the reduced formulation consists of a set of ordinary differential equations describing the wall-normal evolution of the wall-model velocity, the mixture enthalpy and the species mass fractions . The momentum equation reduces to the statement that the wall-normal derivative of the total shear stress, given by the sum of molecular and turbulent contributions, is zero, implying constant total shear stress across the modeled layer. The energy equation accounts for the balance between the work of viscous stresses and the combined molecular and turbulent heat fluxes. The species equations describe the balance between molecular and turbulent diffusion of each species and the corresponding chemical source terms. The mixture thermodynamic properties are evaluated from the local species composition and temperature. In particular, the mixture enthalpy and specific heat at constant pressure are obtained as mass-fraction-weighted sums of the corresponding species properties. Turbulent heat and mass diffusivities are modeled through eddy-diffusivity relations, where the turbulent thermal conductivity is proportional to the turbulent viscosity and inversely proportional to the turbulent Prandtl number, and the turbulent mass diffusivity is proportional to the turbulent viscosity and inversely proportional to the turbulent Schmidt number. The turbulent viscosity μt is obtained using the mixing length model corrected with the Van Driest damping function, while the turbulent Prandtl and Schmidt numbers will be tested either as constants (Prt = 0.9, Sct = 0.7) or through a formulation based on the generalized Reynolds analogy [4]. In the a priori analysis, the DNS database is used to directly extract near-wall gradients, total shear stress distributions and thermodynamic property variations. These DNS-resolved quantities are introduced into the reduced one-dimensional equilibrium wall-model formulation in order to evaluate its assumptions under both perfect-gas and reacting conditions. In this context, the validity of the constant total-stress hypothesis and equilibrium turbulence closure is tested by comparing the wall-model prediction against the exact DNS wall shear stress. In contrast, the a posteriori analysis consists of embedding the wall-model formulation within the simulation framework and evaluating its performance in a fully coupled numerical environment, without directly imposing DNS values. This allows assessment of the model behavior under realistic WMLES-like conditions, whereby the WALE model is employed to compute the subgrid-scale viscosity. The influence of the wall-model exchange location is assessed by varying the matching position within 10 ≤ y+ ≤ 50. For the reacting multi-species Mach 3 cases, turbulence statistics, species mass-fraction distributions and dissociation degree are analyzed. The dissociation degree exhibits a quadratic variation across the channel core, consistent with recent high-enthalpy DNS studies. Compared to perfect-gas conditions, thermochemical non-equilibrium modifies viscosity distributions, temperature gradients and near-wall turbulence structure, directly affecting wall-model assumptions.

5. Conclusions A consistent DNS-based framework for hypersonic wall-model development has been presented. A perfect-gas database of compressible turbulent channel flows has been generated at Mach 1.5 and Mach 3 and both configurations have shown excellent agreement with available reference DNS data, providing verification of the numerical methodology. The extension of the database to multi-species reacting Mach 3 cases is ongoing and will allow a systematic evaluation of thermochemical non-equilibrium effects on near-wall modeling. Based on this DNS reference framework, systematic a priori and a posteriori studies of wall-model assumptions have been performed. The a priori assessment quantifies the validity of equilibrium and constant-stress hypotheses using DNS-resolved near-wall fields, while the a posteriori assessment evaluates the same wall-model formulation when embedded in a fully coupled numerical environment. The present work establishes a coherent link between DNS database generation and physics-based wall-model assessment for future hypersonic Wall-Modeled Large-Eddy Simulations.

References [1] F. Spina, A. J. Smits, and S. K. Robinson. The physics of supersonic turbulent boundary layers. Annual Review of Fluid Mechanics, 1994. https://doi.org/10.1146/annurev.fl.26.010194.001443.

[2] M. Cogo, F. Salvadore, F. Picano, and M. Bernardini. Direct numerical simulation of supersonic and hypersonic turbulent boundary layers at moderate-high reynolds numbers and isothermal wall condition. Journal of Fluid Mechanics, 2022. https://doi.org/10.1017/jfm.2022.574.

[3] D. Passiatore, L. Sciacovelli, P. Cinnella, and G. Pascazio. Thermochemical non-equilibrium effects in turbulent hypersonic boundary layers. Journal of Fluid Mechanics, 2022. https://doi.org/10.1017/jfm.2022.283.

[4] H. Huang, H. Su, Q. Guo, P. Liu, X.Yuan, and J. Cai. Ordinary-differential-equation-based wall model for large eddy simulation of hypersonic high enthalpy turbulent boundary layer. Physics of Fluids, 2025. https://doi.org/10.1063/5.0295932.

[5] G. Pascazio, D. Ninni, F. Bonelli, and G. Colonna. Hypersonic flows with detailed state-to-state kinetics using a GPU cluster. In Plasma Modeling (Second Edition): Methods and applications. 2022. https://doi.org/10.1088/978-0-7503-3559-1ch10.

[6] D. Modesti and S. Pirozzoli. Reynolds and mach number effects in compressible turbulent channel flow. International Journal of Heat and Fluid Flow, 2016. https://doi.org/10.1016/j.ijheatfluidflow.2016.01.007.

ABSTRACT. We analyze the excitation of second Mack-mode disturbances in the boundary layer over the HIFLIER flight vehicle during its ascent. The various sources of disturbances considered relevant to free flight, including free stream turbulence and particulates, are analyzed within the CHAMPS CFD framework. A combination of high-fidelity simulation capabilities and reduced-order modeling are employed to determine which sources were most likely responsible for initiating transition to turbulence over the straight cone flying at approximately Mach 6.

ABSTRACT. High-order numerical methods offer greater accuracy and efficiency than traditional low-order approaches, but their practical adoption can be limited by stability and meshing challenges. Ensuring nonlinear stability is key to making these methods robust and reliable for real-world engineering problems, especially in situations involving moving meshes—such as fluid-structure interactions, moving parts, and adaptive grids. The compressible Navier–Stokes (CNS) equations possess a convex entropy function, related to thermodynamic entropy. This function acts as an admissibility criterion for identifying physically relevant weak solutions and provides a means of demonstrating nonlinear (entropy) stability by bounding in the L2 norm the solution. Maintaining entropy stability for CNS equations has been a major focus of research. Recent advances show entropy stability can be preserved even as computational grids deform, using a discontinuous Galerkin split form of the arbitrary Lagrangian Eulerian (ALE) approach. This work extends these results, proving entropy stability for moving wall boundary conditions and demonstrating the effectiveness of high-order entropy-stable DG methods in large-scale simulations with mesh motion.

ABSTRACT. Numerical discretization of spatial derivatives in partial differential equations inherently introduces errors due to the replacement of the continuous model by its discrete counterpart. At the operator level, discrete derivative approximations can be interpreted as exact derivatives acting on implicitly filtered representations of the solution. Based on this concept of numerically induced filtering, the Approximate Deconvolution Discretization (ADD) method was previously proposed to compensate for discretization-induced filtering effects. However, the original ADD formulation was restricted to periodic boundary conditions.

In this work, we extend the ADD framework to problems with non-periodic boundary conditions. The induced filter is constructed locally using Lagrange interpolation on asymmetric stencils near the boundaries, while maintaining the formal order of accuracy. The ADD procedure is applied to the first-derivative operator, whereas the second derivative is discretized using a stable high-order compact scheme in order to isolate the stability properties of the modified first-derivative operator.

The proposed formulation is assessed using the two-dimensional vorticity–stream function formulation of the Navier–Stokes equations and the regularized shear-driven cavity flow (Burggraf flow), for which an analytical solution is available. Stability is analyzed through eigenvalue spectra of the modified operators. Although direct inversion of the Lagrange-based induced filter leads to operators that are theoretically unstable for certain stencil sizes, numerical simulations demonstrate that stable and highly accurate solutions can still be obtained for low and moderate Reynolds numbers, and for higher Reynolds numbers on sufficiently fine grids. In these cases, viscous dissipation suppresses the growth of unstable numerical modes, leading to conditionally stable behavior despite operator-level instability.

The results show significant accuracy improvements compared to standard second-order finite difference discretization. Ongoing work focuses on identifying stability limits of the approach and on developing strategies to mitigate near-boundary instabilities.

ABSTRACT. Multicomponent real-fluid flows are fundamental to a wide range of engineering and natural systems, such as high-speed vehicles, energy conversion devices, chemical and process engineering, and extreme weather processes. Accurate representation of real fluid flows is critical for advancing predictive simulation tools that support safer, more efficient, and more sustainable engineering designs. This work develops three-dimensional (3D) high-order numerical methods for multicomponent flow simulation with a real-fluid equation of state (EOS) representation. The diffuse interface concept is adopted to approximate flow discontinuities, such as shock waves and material interfaces. As a result, the localized artificial viscosity method is well suited for flow discontinuity capturing as it can spread a discontinuity over a finite thickness, which is controlled by the diffusivity strength. A conservative thermodynamic calculation approach is used to evaluate thermal states of real fluids, which can effectively eliminate the so-called spurious pressure oscillation, even for real fluids with highly nonlinear EOS, although it has to increase the computational cost spent on thermal state evaluation. The real-fluid effects on the Richtmyer–Meshkov instability (RMI) are then studied with high-order numerics. It finds that high-order schemes significantly outperform their low-order counterparts.

ABSTRACT. We propose a high-order discontinuous Galerkin (DG) method for solving the eikonal equation with improved stability and efficiency. Conventional high-order solvers rely on nonlinear equations with auxiliary variables, increasing complexity and cost. We extend a linearized vanishing-viscosity algorithm to a high-order DG method, but direct extension leads to projection-induced error amplification and oscillations near singularities. To address this, we terminate vanishing-viscosity iterations at an optimal stage and switch to an inviscid formulation to prevent error amplification in smooth regions. A robust criterion is also proposed to detect singular regions generated by characteristic front collisions. The Singular Gradient Regularization Method (SGRM) then applies viscosity only in singular regions while keeping smooth regions inviscid. A region-decoupled flux prevents viscosity contamination of smooth regions. The method uses a scalar-only formulation, avoiding auxiliary gradient equations and reducing computational cost. Extensive numerical experiments on complex 2D and 3D geometries demonstrate high-order accuracy and stable behavior near singular gradients. Integration into high-order RANS and hybrid RANS–ILES solvers confirms its effectiveness for practical wall distance computations.

ABSTRACT. Urban areas currently support more than 55% of the global population, a figure projected to reach 68% by 2050. High population density and restricted evacuation options elevate the risk of human exposure to airborne contaminants. Consequently, developing rapid response capabilities for the accidental or intentional release of chemical, biological, or radiological substances is a critical challenge. In these scenarios, the application of physics-based urban wind and dispersion models is essential for emergency management protocols and decision-making frameworks.

Ensuring the rapid turnaround required for emergency response necessitates models that are both highly automated and computationally efficient. Cartesian grid methods with embedded boundaries address the automation requirement by facilitating rapid, fully automated mesh generation for complex geometries. Beyond streamlined meshing, these methods offer flow-based adaptation, straightforward numerical formulations, and low numerical dissipation. Furthermore, the inherent regularity of Cartesian grids enables efficient vectorisation and is well-suited for high-performance parallel architectures such as GPUs. To complement these advantages and ensure rapid convergence, multi-level methods are employed to mitigate the deteriorating convergence rates typically exhibited by iterative solvers, such as the SIMPLE algorithm, as grid density increases. Specifically, the non-linear Full Approximation Scheme (FAS) accelerates convergence by efficiently resolving low-frequency error components on coarser grids. Previous investigations have demonstrated the effectiveness of the FAS multigrid approach for a variety of complex problems, including urban windflows.

To meet these demands, we developed a Cartesian grid solver that combines a cut-cell embedded boundary method with a multilevel algorithm to achieve both efficiency and automation. Developed using the AMReX framework, the method solves the steady-state, incompressible Reynolds-Averaged Navier-Stokes (RANS) equations using the SIMPLE algorithm. Turbulence is modelled with the Wilcox k-omega (2006) model, utilising a high-Reynolds number wall treatment for cut-cells. To automate the mesh generation, we use flow based mesh adaptation using sensors to detect areas that need refinement. In this short abstract, we outline briefly the implementation of the high-Re wall treatment and present some validation results.

ABSTRACT. The mechanical interactions between ultrasound-stimulated microbubbles (USMB) and cell membranes are a critical phenomenon in targeted drug delivery and immunotherapeutic strategies. The acoustically driven oscillation of the microbubbles has been shown to induce sonoporation, a transient increase in cell membrane porosity that enhances site-specific drug delivery. In contrast, inertial cavitation can irreversibly damage and kill cells, which should be avoided during USMB-based therapeutic and diagnostic procedures. Direct observation of ultrasound-stimulated microbubble (USMB) interactions with cells remains challenging due to the fast timescales (microseconds) and small dimensions (micrometers) involved. Nevertheless, optimizing ultrasound-based diagnostic and therapeutic techniques requires a detailed understanding of these interactions at the cellular interface. A comprehensive knowledge of how USMBs behave in proximity to cells is essential for enhancing both their clinical efficacy and safety. To this end, we have developed an interface-resolved axisymmetric computational model to investigate these interactions. In our model, the cell is represented as a droplet of specific viscosity and density suspended within a carrier fluid. A phase-field method is used to model the degradation of the intracellular structure due to the cavitation-induced strain. Our model captures the essential bubble-droplet interactions and interfacial mechanics necessary to characterize USMB-mediated cellular effects. The results show that the stably oscillating microbubbles translate toward the denser and more viscous droplets in the carrier fluid. As the bubble gets closer to the cell, the flow induced by its oscillation periodically deforms the droplet interface. Such deformation of the cell boundaries in cell-bubble interaction may instigate the formation of a temporary pore at the cell membranes, enhancing drug delivery through sonoporation. Furthermore, the acoustically driven bubbles are seen to form high-speed micro jets toward the droplets. The impingement of high-velocity microjets during the inertial cavitation of USMBs can cause irreversible damage, leading to cell lysis—an undesirable outcome in ultrasound-mediated drug delivery. The degradation model will track the damage caused to the intracellular structure due to microjetting. However, jetting collapse can be moderated by adjusting the ultrasound frequency. Our calculations demonstrate that microjet formation is suppressed at higher acoustic frequencies. These results indicate that the current model provides critical insights into cell-bubble interactions, serving as a framework for the development and optimization of USMB-based diagnostic and therapeutic techniques.

ABSTRACT. This work presents advancements in Computational Multi-Fluid Dynamics (CmFD) by proposing a thermally driven Marangoni flow model in a Sharp Interface Level Set Method on a co-located grid (SI-LSMcol). Further, the associated formulation, discretization, and implementation of the interfacial jump condition, essential on a co-located grid, are presented. Finally, validation of the in-house two-phase flow solver is validated for a range of benchmark CmFD problems, including thermocapillary single-fluid Stokes flow in an open cavity, linear shear two-fluid flow, and thermocapillary two-fluid flow in a closed cavity.

ABSTRACT. We present a systematically second-order accurate method for the local evaluation of hydrodynamic loads in immersed boundary simulations for fluid–structure interaction applications. The approach is developed within a finite-difference direct numerical simulation solver for the incompressible Navier–Stokes equations, coupled with a sharp-interface immersed boundary formulation.

The work identifies the discrete mechanism responsible for the well-known degradation of viscous stress accuracy near immersed interfaces. The loss of convergence is traced to a breakdown of local discrete error regularity at boundary-affected grid points, where stencil values no longer share a consistent discretization error structure and second-order cancellations fail.

To prevent this degradation, a probe-based reconstruction strategy is introduced in which pressure and velocity gradients are evaluated exclusively at regular fluid locations and subsequently extrapolated to the immersed surface. Probe placement explicitly avoids boundary-affected points, ensuring consistent second-order accuracy without empirical tuning of probe distances.

The method is verified through manufactured-solution tests and validated against reference numerical and experimental benchmarks, including the flow past a sphere in free stream and sphere sedimentation under gravity. Second-order convergence of reconstructed surface tractions is demonstrated.

ABSTRACT. Rocket-Based Combined Cycle (RBCC) engine integrates a high thrust-to-weight ratio rocket and a high specific impulse dual-mode ramjet within a single flow path, enabling efficient and economical flight across a wide speed range. This makes RBCC a key propulsion solution for future reusable space transportation vehicles. The combustion chamber of RBCC engine is persistently exposed to extreme conditions involving high temperature, high heat flux, and high-speed flow. An efficient Thermal Protection System (TPS) is essential to maintain structural temperatures within material limits, and thermal protection technology has become a critical bottleneck in engine development. Among various active cooling strategies, regenerative cooling has attracted significant research and engineering interest due to its high efficiency, stable wall temperature control, improved energy utilization, and potential to enhance fuel combustion performance, thereby supporting long-duration and reusable engine operation. Conventional regenerative cooling systems typically employ straight rectangular channels, favored for their structural simplicity and low cost. However, under multi-mode operating conditions, strong coupling between supersonic combustion and regenerative cooling results in a non-uniform distribution of thermal environment and heat load. Moreover, hydrocarbon fuels often operate under supercritical conditions, where drastic variations in thermophysical properties can induce flow maldistribution and localized overheating, potentially leading to engine failure. Therefore, optimizing cooling channel design is critical to improve flow distribution and overall cooling performance. Topology optimization has become a key methodology for designing channel configurations that achieve optimal performance under specified constraints, with substantial progress in fluid flow and conjugate heat transfer applications[1][2]. Yu et al. [3]employed the moving morphable components (MMC) method to design heat dissipation channels, achieving improved stability compared to density-based methods and reducing the formation of unmanufacturable micro-channels. In open-source platform development, U. Nilsson et al.[4] implemented a topology optimization framework in OpenFOAM using the continuous adjoint method for sensitivity analysis, systematically investigating adjoint equation derivation and discretization of adjoint boundary conditions. Matsumori et al. [5]optimized thermo-fluid channel configurations with the objective of maximizing heat transfer, proposing a design methodology for heat exchanger channels under constant power input and evaluating Reynolds number effects on topological complexity. Research from the Technical University of Denmark (DTU)[6] explored simplified models for fluid topology optimization in natural convection problems, showing that although boundary layer effects are not captured, such models provide effective initial guesses for high-fidelity optimization, reducing computational cost. Existing studies and our previous study[7] indicate that topology-optimized regenerative cooling channels can overcome key limitations of conventional straight channels, including flow maldistribution and localized heat transfer deficiencies. However, state-of -art researches focus on thermal management of cooling structure with small temperature variations. Large thermophysical property changes due to temperature increase from 300K to 900K in practical engine cooling applications make topology optimization of regenerative cooling channels much more difficult. In this paper, topology optimization model considering large property variation for supercritical hydrocarbon fuel via CoolProp is applied to the regenerative cooling channels of a typical ramjet combustor wall. Three-dimensional numerical simulations are conducted to compare flow and heat transfer characteristics between the optimized and conventional straight channels, validating that topology optimization enables the design of lighter regenerative cooling structures with enhanced heat dissipation efficiency. Figure 1 shows a schematic diagram of a typical structural component of a ramjet combustor. Its wall structure is extracted and simplified into the topology optimization domain for the present study. Owing to the symmetry of the physical problem, half of the full topology optimization domain is adopted as the actual optimization domain. The design domain consists of a 100 mm × 13.75 mm rectangular region and six inlets and outlets with a dimension of 4 mm × 2 mm, where the spacing between adjacent inlets and outlets is 2 mm. Each inlet and outlet is located 1 mm away from the symmetry plane and 2.75 mm away from the top surface of the design domain. A volumetric heat source with a magnitude of 1×10⁸ W/m³ is applied within the design domain. The inlet velocity is set to 0.05 m/s, and the inlet temperature is set to 300K. For topology optimization problems involving thermo-fluid coupling, the density method is commonly adopted for optimization. A uniform volumetric heat source is applied within the design domain, whose magnitude is determined by the average heat flux of the heated wall. With prescribed energy dissipation constraint and volume constraint, the optimal spatial distribution of fluid-solid material elements is solved via numerical calculation, so as to minimize the objective function of the average temperature in the domain. In the present work, the volume constraint is prescribed at 0.6. Based on the density method, the topology optimization formulation for the thermo-fluid problem with the objective of minimizing the average temperature can be expressed as follows: In the present work, the high-temperature alloy GH3128 and n-decane are selected as the solid and fluid materials, respectively. Given the drastic variations of n-decane’s thermophysical properties with pressure and temperature in the channel during regenerative cooling, accurate calculation of its thermophysical properties is required, for which polynomial fitting and tabulation with interpolation are the two common approaches. Herein, the tabulation with interpolation method is employed. The n-decane thermophysical property table is established via CoolProp at 3 K temperature intervals and 1 kPa pressure intervals without considering fuel pyrolysis, and the fluid thermophysical parameters in the numerical simulation are acquired by linear interpolation. The final topology optimization result from the calculation is shown in Fig. 2. The ramjet combustor wall is typically divided into an inner wall, a channel layer, and an outer wall. The channel layer is formed by extruding the topology optimization result, and both the inner and outer walls are fabricated from the high-temperature alloy GH3128, to establish the three-dimensional combustor wall model. Numerical simulations of flow and heat transfer are conducted on the three-dimensional topology-optimized regenerative cooling structure. The standard k-ε model is selected for the Reynolds-Averaged Navier-Stokes (RANS) calculations, which has been widely validated to accurately predict the flow and heat transfer of supercritical-pressure hydrocarbon fuel in horizontal channels. The Coupled algorithm is used for pressure-velocity coupling, with the convection and diffusion terms discretized by the second-order upwind and second-order central differencing schemes, respectively. Given the complex species and reaction processes involved in the endothermic pyrolysis of kerosene (a multi-component hydrocarbon fuel mixture), a surrogate component model is adopted for numerical calculation. A simplified three-step global reaction mechanism with 17 species is employed to simulate the endothermic pyrolysis of kerosene under supercritical pressures, with C₁₁.₈₅H₂₃.₈₂ as the kerosene surrogate. The simulation results of three-dimensional flow and heat transfer considering pyrolysis are shown in Fig. 3. The design of engine regenerative cooling channel via the topology optimization method holds great promise for improving flow distribution uniformity and heat transfer efficiency, thus enabling more efficient wall thermal protection. In the present work, three-dimensional numerical simulations are conducted to evaluate the cooling performance of the topology-optimized cooling panel under the operating heat flux conditions of an RBCC engine, with the main conclusions drawn as follows: (1) In terms of structural weight, the solid domain volume fraction of the channel layer is 43% for the topology-optimized configuration, compared to 58.4% for the conventional straight configuration, corresponding to a 26.3% reduction in channel layer weight. (2) Regarding thermal characteristics, the topology-optimized configuration exhibits an average temperature rise of 698 K on the heated wall, whereas that of the conventional straight configuration is 992 K, achieving a 29.6% reduction. The maximum temperature rise on the heated wall reaches 799 K for the topology-optimized configuration, in contrast to 1099 K for the conventional straight configuration, representing a 27.3% decrease. (3) Under identical fluid volume fraction conditions, the topology-optimized channel achieves a 28.1% reduction in average temperature rise and a 25.2% reduction in maximum temperature rise on the heated wall compared to the conventional straight channel. This suggests that the reduction in hot wall temperature achieved by the topology-optimized channel is primarily attributable to the channel topology optimization itself, rather than an increase in fluid volume fraction. (4) The reduction in wall temperature of the optimized regenerative cooling channel becomes more pronounced under low inlet temperature conditions. The present study demonstrates that topology optimization serves as a feasible and effective computational design tool for regenerative cooling systems in ramjet engines. Future efforts will be directed toward the additive manufacturing of the topology-optimized cooling structure, followed by ground-based direct-connect thermal qualification tests of the ramjet combustor, with the aim of further facilitating the engineering implementation of topology optimization technology in the aerospace propulsion field.

ABSTRACT. Adjoint-based shape optimization in turbulent flows is developed by explicitly treating turbulent fluxes as forcing terms in the governing equations. Conventional RANS-based adjoint methods rely on eddy-viscosity models, which assume isotropic turbulence and may limit accuracy, particularly near walls. To overcome this limitation, turbulent statistics are first obtained from direct numerical simulation, and the Reynolds stresses and turbulent heat fluxes are incorporated into the adjoint framework as virtual forcing terms under the frozen turbulence hypothesis. The method is applied to the multi-objective optimization of a pin-fin configuration in turbulent channel flow, aiming to reduce drag while enhancing heat transfer. Sensitivity information is derived and used to update the geometry through a level-set method. The optimized design achieves a 13% reduction in drag and a 4% increase in heat transfer, demonstrating the effectiveness of the proposed framework for accurate turbulent shape optimization.

ABSTRACT. 1. Introduction Modern computational fluid dynamics (CFD) workflows increasingly rely on large numbers of parameterized simulation runs to support design exploration, uncertainty quantification (UQ), Non-deterministic Evaluation (NDE) and verification activities. Typical examples include angle-of-attack sweeps, grid refinement studies, and repeated analyses under varying operating conditions. While high-performance computing (HPC) platforms enable these studies to be executed efficiently, ensuring that results remain reproducible, auditable, and reusable over time remains a persistent challenge. A well-defined workflow for one instance of the parameterized study, that may even be tracked through version control, addresses how information is derived through a process since it essentially describes the utilized algorithms. Provenance addresses what information is derived through a process since it describes relationships within a set of files for each instance of the parameterized study including reference to the utilized workflow. Provenance, rather than workflow alone, enables an engineer to trace to specific files to investigate interesting nuances to gain insight into the overall results as well as investigate questionable results that drive CFD model refinement (e.g. flow visualization from volume solution files corresponding to one point on a UQ curve).

The current study addresses the problem of CFD results traceability and reproducibility by treating the set of files utilized in a process and associated relationships between them as a first-class computational object. The management, or data control, of a process’ files occurs through a provenance-aware execution method that explicitly records the relationships among simulation parameters, solver runs, generated artifacts, and post-processing steps as workflows execute. Furthermore, file access is virtualized, simplifying review of the many file artifacts used and created during workflow execution. The approach is demonstrated using a parametric CFD workflow based on the XFOIL panel method, representative of the repeated analyses commonly performed in preliminary aerodynamic design and UQ studies. The value of the provenance-aware method can be clearly understood for efficiently working through the use cases of large parameterized studies outlined below.

1. Traceability: Trace specific process instance files to: a. Find source file(s) for a runtime or logic error b. Understand how a CFD model should be refined c. Investigate specific run instances corresponding to interesting nuances of the overall results 2. Reproducibility: Invert retrospective provenance (what actually occurred) into a precise and detailed reproduction workflow to regenerate results rather than hoping for equivalent execution of a pre-defined workflow, which is prospective provenance (what was planned to occur)

2. Methodology 2.1 Workflow-Centric CFD Execution Model In the proposed approach, a CFD study is represented as an executable workflow composed of well-defined stages and run units, which defines how information will be derived. Each run corresponds to a solver invocation with a specific set of input parameters, while stages group related runs and post-processing steps. Key workflow elements include airfoil geometries, run parameters (e.g., Mach number, angle of attack), solver input and output files, scripts, and derived results such as plots or aggregated quantities of interest. By expressing these elements explicitly, the workflow structure becomes independent of any particular file system layout or scripting convention. 2.2 Provenance Capture During Execution Provenance information is captured automatically as workflows execute, recording both the intended workflow structure and the actual execution history, which defines what information is derived. This includes the generation and consumption of data artifacts, the parameters used for each solver run, the software components involved, and the ordering of workflow activities. Capturing provenance at execution time ensures that the full chain of derivation for any result can be reconstructed without relying on external documentation or manual bookkeeping. The capture process is designed to be open and non-intrusive, allowing engineers to run analyses using familiar tools and HPC environments while provenance is recorded transparently. 2.3 XFOIL-Based Parametric and UQ Workflow The methodology is demonstrated using a panel-method aerodynamic analysis implemented with MIT’s XFOIL code. In the example workflow, the airfoil geometry and Mach number are held fixed while the angle of attack is varied over a prescribed range, producing a sequence of solver runs whose results are subsequently consolidated and post-processed to generate aerodynamic polars. The second example is a UQ workflow, in which XFOIL runs are used to create a surrogate model for statistical sampling. These activities are used to perform sensitivity analysis (SOBEL indices) and compute confidence metrics. These examples are representative of common parametric sweeps and UQ studies in CFD practice, where tens or hundreds of solver runs must be managed consistently. 2.4 Results and Reproducibility Analysis Captured provenance allows the content for any generated artifact to be traced in terms of how and what information was derived or reproduced based on the exact parameters, solver version, and intermediate results involved in its creation. The provenance representation supports direct re-execution of the entire analysis suite or selected subsets with modified parameters, enabling controlled comparisons between runs. Differences between workflow execution instances—such as parameter changes or altered solver settings—can be identified explicitly, providing a concrete basis for in-depth investigations, reproducibility assessment and auditability in CFD studies. Organizing artifacts via their workflow provenance provides an additional benefit as users can access file artifacts directly by clicking on thumbnails in the provenance timeline. This eliminates the distraction of tracking down file names and folder paths when seeking details within complex workflows. 3. Conclusions This work demonstrates that explicit provenance capture integrated with workflow execution provides a practical foundation for efficient in-depth investigations to better understand a broad set of results and reproducible CFD analysis, particularly for parametric and uncertainty-driven studies. No additional burden is on the analyst for reliable rerunning, comparison, and auditing of simulation results when using this approach since the total provenance relationship structure spans the pre-defined workflow, all artifacts recorded during execution, users, computers and software, and execution activities that occur over time. The XFOIL-based examples will illustrate how common CFD workflows can be made reproducible and transparent, with direct relevance to verification, benchmarking, and UQ applications. Future work will extend this approach to more complex CFD solvers and multi-physics workflows.

ABSTRACT. 1. Introduction

The design and development of low-boom supersonic aircraft for overland flight has attracted significant interest from both public and private sectors. NASA’s Quesst mission, for example, aims to demonstrate the feasibility of quiet supersonic travel using the X-59 low-boom flight demonstrator, where public response to sonic boom loudness is a critical metric. As the perceived loudness of a sonic boom at ground level is strongly influenced by aircraft shape [1], high-fidelity computational fluid dynamics (CFD)-based design optimization is a well-positioned tool for next-generation supersonic vehicle design.

Gradient-based optimization provides a tractable approach for exploring design spaces efficiently. Early work by Rallabhandi et al. [2] introduced adjoint-based optimization for sonic boom minimization using the CFD solver FUN3D coupled with the noise propagation tool sBOOM [3], leveraging the Euler equations to model the aerodynamics in the nearfield. Subsequently, Rodriguez et al. [4] integrated the Cartesian-based Euler solver Cart3D with sBOOM for output-based mesh adaptation and noise-based shape optimization. Building on these foundations, this work applies both Euler and Reynolds-averaged Navier-Stokes (RANS)-based analysis on curvilinear overset grids for nearfield aerodynamics, coupled with sBOOM for farfield noise propagation, within a multidisciplinary optimization framework. This approach enables accurate modeling of viscous effects while maintaining flexibility for future extensions to additional disciplines.

2. Methodology

2.1 Flow and Sonic Boom Propagation Solvers

The ground-level sonic boom predictions performed in this work rely on a two-step process. Nearfield aerodynamics are modeled using the Launch, Ascent, and Vehicle Aerodynamics (LAVA) Curvilinear solver [5], which can solve the Euler and RANS equations on structured curvilinear overset grids. For RANS-based analyses, the negative variant of the Spalart-Allmaras (SA) turbulence model is applied. A discrete adjoint solver is also applied for sensitivity analysis. The input waveform for farfield propagation is constructed by interpolating nearfield overpressure values onto a propagation path.

Farfield sonic boom propagation is then performed using sBOOM [3], which solves the augmented Burgers equation accounting for nonlinearity and atmospheric effects. sBOOM also computes noise metrics and provides adjoint capabilities for sensitivity analysis with respect to the input waveform.

2.2 Optimization Problem and Gradient Calculations

The design optimization problems considered in this work seek to minimize ground-level sonic boom metrics by controlling geometric design variables. Due to the partitioned approach adopted for modeling ground-level noise, where nearfield aerodynamics and farfield propagation are solved with separate methods, a coupled adjoint methodology is required to compute sensitivities efficiently. The coupled adjoint formulation used here follows the methodology introduced by Rallabhandi et al. [2].

2.3 Geometry Parameterization and Grid Deformation

Geometries are parameterized using either OpenVSP or free-form deformation (FFD) volumes through the open-source pyGeo package [6], enabling control of user-specified design variables for shape changes during optimization.

For nearfield aerodynamics, we employ structured curvilinear overset grids comprising a near-body grid that conforms to the surface and an off-body grid that extends to the nearfield domain boundaries. In supersonic applications, the off-body grid is constructed to align with the freestream Mach angle, which reduces numerical dissipation along the dominant flow characteristics. To preserve this benefit across design iterations, only the near-body grid is deformed in response to geometric changes, while the off-body grid remains fixed and aligned with the Mach angle. This is demonstrated in Figure 1.

At each design iteration, overset connectivity and interpolation weights are recomputed to reflect the updated near-body topology and any blanking/unblanking changes in the off-body grid. Resulting changes in overset connectivity and interpolation are accounted for within gradient calculations, so that gradients with respect to design variables remain accurate under near-body deformations with fixed off-body grids.

3. Results

We present two results to demonstrate the applicability of design optimization for sonic boom ground noise minimization. The first case is the optimization of a two-dimensional diamond airfoil to minimize the mean square of the ground-level pressure signature. The second case is the optimization of the wing twist distribution of the JAXA wing-body configuration for the minimization of the A-weighted sound exposure level (ASEL).

3.1 Diamond Airfoil Optimization

The shape of a diamond airfoil is optimized in order to reduce the mean square of the ground-level pressure signature evaluated directly under the flight path. This simplified two-dimensional case allows us to validate the optimization framework and sensitivity calculations before applying the methodology to more complex three-dimensional configurations. The shape of the diamond is parameterized using OpenVSP such that the inflection point on the lower surface can be adjusted vertically, while the surface segments before and after this point remain linear. The Euler equations are used to model the aerodynamics in nearfield at a Mach number of 2.0 and an angle of attack of 0◦.

The initial and optimized geometries, nearfield pressure perturbations, and ground-level pressure signatures are shown in Figure 2. From these results, we can see that the optimizer almost completely flattens the lower surface of the diamond airfoil. This is done in order to minimize any “turning” in the flow traveling under the airfoil, thus avoiding the introduction of oblique shocks or expansion waves. This results in a nearly zero ground-level pressure signature for the optimized design, as seen in Figure 2c. The initial and optimized mean-square values of the ground-level pressure signatures are 0.117 and 1.01 × 10−7, respectively.

3.2 JAXA Wing-Body Optimization

The JAXA Wing-Body geometry from the 2nd Sonic Boom Prediction Workshop is optimized in order to reduce the ground-level ASEL metric evaluated directly under the flight path. The geometry is parameterized using the FFD volume shown in Figure 3, where the red spheres indicate the FFD control points. In this optimization, the four outboard spanwise stations of the FFD are allowed to rotate about the midchord of the wing, providing control over the wing’s twist distribution. The RANS equations with the SA turbulence model are used to model the aerodynamics in the nearfield, at a Mach number of 1.6, angle of attack of 0◦, and Reynolds number per meter of 5.7 million.

Figure 4 presents a comparison between the initial and optimized geometries as well as the initial and optimized ground-level overpressure signatures. We observe that the optimizer has created some slight oscillations in the twist distribution of the wing, which has led to a change in the aft portion of the ground signature. These changes in signature have reduced the ground-level ASEL value from 62.1 dB for the initial geometry to 59.6 dB for the optimized geometry. Thus, the optimizer has successfully reduced ASEL by 2.5 dB.

4. Conclusions

This work demonstrates a gradient-based optimization framework for minimizing ground-level sonic boom metrics using a partitioned approach that couples RANS-based nearfield analysis with farfield sonic boom propagation. The methodology leverages curvilinear overset grids to maintain accuracy during geometric changes and employs a coupled adjoint strategy to compute sensitivities efficiently across nearfield and farfield domains. Initial results confirm the effectiveness of the approach and its potential for application to supersonic aircraft designs.

The final paper will contain shape optimization for the C25D configuration from the 2nd Sonic Boom Prediction Workshop, a more complex geometry than the JAXA Wing-Body presented herein.

References

[1] Joseph Pawlowski, David Graham, Charles Boccadoro, Peter Coen, and Domenic Maglieri. Origins and overview of the shaped sonic boom demonstration program. In 43rd AIAA Aerospace Sciences Meeting and Exhibit. AIAA, January 2005. AIAA 2005-5.

[2] Sriram K. Rallabhandi, Eric J. Nielsen, and Boris Diskin. Sonic-boom mitigation through aircraft design and adjoint methodology. Journal of Aircraft, 51(2):502–510, March 2014.

[3] Sriram K. Rallabhandi, Marian Nemec, and Michael J. Aftosmis. Recent enhancements to modeling sonic boom propagation using augmented burgers’ equation. In AIAA AVIATION 2023 Forum. AIAA, June 2023. AIAA 2023-3727.

[4] David L. Rodriguez, Michael J. Aftosmis, Marian Nemec, and Sriram Rallabhandi. Sonic boom ground noise minimization via the adjoint method. In AIAA SCITECH 2025 Forum. AIAA, January 2025.

[5] Cetin C. Kiris, Jeffrey A. Housman, Michael F. Barad, Christoph Brehm, Emre Sozer, and Shayan Moini-Yekta. Computational framework for launch, ascent, and vehicle aerodynamics (LAVA). Aerospace Science and Technology, 55:189–219, August 2016.

[6] Hannah M. Hajdik, Anil Yildirim, Neil Wu, Benjamin J. Brelje, Sabet Seraj, Marco Mangano, Joshua L. Anibal, Eirikur Jonsson, Eytan J. Adler, Charles A. Mader, Gaetan K. W. Kenway, and Joaquim R. R. A. Martins. pyGeo: A geometry package for multidisciplinary design optimization. Journal of Open Source Software, 8(87), July 2023.